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ABSTRACT 


The Third Quarterly Technical Report for the LASA Experimental Signal 
Processing System identifies the effort expended to provide the hardware and 
software in support of research and development directed toward the study of 
seismic signal processing. It also delineates work tasks planned for the next 
quarter. This document presents detailed information related to machine con¬ 
figurations, time delay correlation, event location accuracy, optimum process¬ 
ing, fast Fourier transform, array design, steering delay library, magnitude 
estimation, and travel time characterization. 
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Section 1 


INTRODUCTION 


The work reported herein was performed under the "LASA Experimental 
Signal Processing System n Contract Number F19628-67-C-0198, and is a con¬ 
tinuation find extension of the "Large Aperture Seismic Array Signal Processing 
Study,"1 Contract Number SD~296,and the n Large Aperture Seismic Array 
Signal Processing Communications and Simulation Study, "2,3,4 Contract Number 
A FI 9(62 8)-5948. 

The purpose of this contract is to design, develop, and implement the LASA 
Experimental Signal Processing System (ESPS), including the hardware and soft¬ 
ware, to provide an experimental capability to: 

a. Evaluate performance of the system in accordance with system 
requirements 

b. Demonstrate the capability to meet basic LASA signal processing 
objectives 

c. Conduct research to develop means of improving and extending 
the capability of the system 

d. Perform seismic and signal processing experiments of interest. 

5 0 

During the first' and second quarters, effort was primarily directed 
toward establishing the Seismic Array Analysis Center (SAAC), Washington, D.C., 
installation and operation of the Detection and Event Processors as System/360 
Model G and II general purpose digital computers, respectively, and the con¬ 
tinuation of system studies in seismic signal processing. 

Work in the present quarter has focused on such signal processing studies 
as time delay correlation, event location accuracy, optimum processing, fixed- 
point Fourier transform, and array design programming for system development. 
In addition, activity in event magnitude estimating, time delay library formula¬ 
tion, and travel time characterization supporting event processing has been 
initiated. 
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Seetion 2 


RESULTS AND SUMMARY OF WORK ACCOMPLISHMENTS 


The following subsections correspond to the functional disciplines of opera¬ 
tions, systems, and programming. 


2.1 OPERATIONS 

The principal effort during this quarter was expended on the maehine sys¬ 
tem, Experimental Display and computer operations. 

Installation and checkout of the SAAC hardware configuration was completed. 
The current and planned configurations are described in Appendix I. The IBM 
2250-1 CRT Display was installed, and debug of the software support packages 
was initiated. 

Modification and integration of the Experimental Display was completed for 
both the off-line and on-line modes. The strip chart recorder interface, 
reported previously has been deferred due to priority work in other areas. 


2.2 SYSTEMS 

During this quarter, effort pertinent to overall process improvement and 
event processing structure has been undertaken. 

As reported in Appendix II-1, a method of automatically determining time 
delays from the event data has been formulated. Preliminary testing indicates 
process stability and that the resulting set of time delays is consistent even when 
the signal-to-noise ratio at the seismometer level is near unity. 

The effects of time delay and system calibration errors on large array 
event location aecuraey are presented in Appendix II-2. It is shown that although 
the processed output record signal-to-noise ratio and system calibration errors 
contribute to location error, the significant location errors are due to array 
aperture and incorrect time delays. 
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A study to determine the effectiveness of fidelity optimum processing to 
cJiscriminate against one of two simultaneously occurring events has been under¬ 
taken. The results presented in Appendix II-3 relate the effect of signal-to-noise 
ratio and array geometry on processing performance. 

Current interest in the Fast Fourier Transform algorithm has prompted 
investigation of its fixed arithmetic execution precision to determine the appli¬ 
cability of microprogramming. It is shown in Appendix II-4 that the errors in 
fixed point execution can be reasonably bounded by scaling the process when an 
overflow occurs and that microcoded execution speed factors are about ten over 
standard assembly language programming. 

A method of array geometry optimization has been established and tested 
and is described in Appendix II-5. Since the technique employs interelement 
distance dependent correlation coefficient data as input, its principal advan¬ 
tage appears to be comparative in nature. 

In the phase delay activity, effort was concentrated on defining process 
structure, automatic acquisition, anomaly characterization, and subarray beam 
assignment. 

With respect to process structure and automatic acquisition, the chief 
accomplishment has been the development of a library organization for handling 
the phase delays for both detection and event processing. Progress on this 
effort is described in Appendix III-1. 

A procedure based on the use of sets of correction factors where each phase 
delay set is valid over a specified geographic area called a region has been 
developed. The procedure also involves the use of interpolation techniques to 
calculate corrections relative to areas between regions. Appendix III—2 describes 
the regions and their associated correction factors as well as a test which was 
performed to evaluate the effectiveness of the region groupings. The results of 
the test indicate that regional corrections improve system performance. 

The subarray beam assignment investigation addresses the method of cover¬ 
ing an arbitrary location in inverse velocity space by a judicious choice of sub¬ 
array beams. During this quarter, formulas were developed for calculating sub¬ 
array steering losses relative to a specified configuration of preformed beams 
and the associated LASA steering loss resulting from the use of such preformed 
beams. In addition, a rule was specified for computing the beam from each sub¬ 
array which is to be used in forming the LASA beam, when given a desired LASA 
beam location. Details regarding the subarray assignment investigation may be 
found in Appendix III-3. 

In Appendix IV-1, the empirical relationships between the classical event 
magnitude measurement and the kinetic energy "density 1 ' method of magnitude 
estimation are given. 
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A significant portion of event processing assumes the use of travel time 
data to sort wave arrivals into consistent arrival families and to use this data 
in estimating event location. The travel time information in its original tabular 
form is bulky and would require interrelation with epicentral range and haz¬ 
ardous differentiation of discrete data. In Appendix IV-2, orthogonal polynomial 
fitting in the least squares error sense was applied to a portion of the Jeffreys- 
Bullen Seismological Tables to explore the existence of explicit functions relating 
body wave travel time and horizontal inverse phase velocity to epicentral event 
range at a known depth. 


2.3 PROGRAMMING 

Programming effort during this quarter was concentrated in Detection Pro¬ 
cessing, Supervisor/Monitor, Plotting Routines, Display Programs, and Experi¬ 
mental Programs. Accomplishments in each of these areas will be presented in 
the following subsections. 


2.3.1 Detection Processing 

At the completion of this quarter's work, three operational versions of the 
Detection Controller were in existence. Two of these were oriented towards a 
real-time operating environment, while the third was tailored to support analysis 
activity. 

The difference between the two real-time versions lies in the area of input/ 
output control. One operates using the standard IBM System/.3(30 Disk Oper¬ 
ating System (DOS), while the other operates using a specially written input/ 
output controller. Some experiments were made to determine if there was any 
timing advantage of one method or the other. The latter was faster by about 
0.5 percent and appears to have no significant effect on the total number of array 
beams that can be performed in a given time interval. 

The data analysis version of the Detection Controller uses standard DOS 
input/output, but it is intended for operation under slightly different conditions. 

In particular, the signal rectification and integration process has been set at 
5 Hz, rather than 1.67 Hz, and the threshold detection and noise integration 
functions have been removed. The capability to produce intermediate tapes of 
subarray beams (unfiltered and filtered), array beams, and short time average 
values has been included in this version. The primary use of this program, 
thus far, has been to produce input tapes for both beam pattern correlation and 
experimental display work. 

All three of these programs have the same basic structure and design 
philosophy. All accept the standard SAAC Edit Format, as input, although they 
do not presently interpret the array status data. A unique feature of these pro¬ 
grams is a storage allocation routine which allows maximum use of available 
core storage for each run of the Detection Controller. This has enabled us to 
process as many as 1024 LASA beams in one run of the program. 
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An additional pi'Ogram capability includes the dispersed subarray beam 
capability (see Appendix III-3). An input card supplied with the steering 
delays for each array beam indicates to the Detection Controller which of the 
subarray beams are to be used for array beamforming. On the basis of this 
information, the delays are properly indexed for use by the beamforming 
microprogram. 


2.3.2 Supervisor/Monitor 

Systems generation and maintenance activity was accented upon availability 
of the System/360 2040H early in the quarter. A multiprogramming Disk Opera¬ 
ting System was provided. Various versions were generated as configuration 
changes took place. At this point, the system contained no alterations, but it was 
generated with a greatly expended supervisor area in anticipation of the future 
core memory requirement. Special attention was given to output spooling tech¬ 
niques which would allow reasonably efficient program development using a 
multisystem including just one 'line printer and one card punch. 


2.3.2.1 Detection Subsystem 

A test which measured the overhead incurred by using a standard Disk 
Operating System supervisor in a detection processing function was conducted. 
The results underscored the advantages in choosing this system as a super¬ 
structure around which the Detection Monitor could be built. 

Work in this quarter also included detail design, coding, and simulated 
testing of the first Detection Monitor package. This package embraced capa¬ 
bilities originally planned as two models. First, the Disk Operating System was 
modified to enable multiprogramming operation in a machine without storage 
protection so that the data acquisition function could operate asynchronously 
with the detection processing function, while having common access to certain 
portions of core memory. The monitor was designed so that inputs could be 
received in real time from any source, with the detection processing programs 
relieved of device-dependent considerations. The same design also enabled the 
setting of mode indicators from external sources to allow dynamic entry by the 
detection processing program into various outputting modes for display purposes. 

The design approach is notable in that the Detection Controller may be inte¬ 
grated with the Detection Monitor into a working system merely by assembling 
macros with the program. The code generated by the macros represents the 
Detection Monitor and becomes a separate program only at the time the Detec¬ 
tion Subsystem is loaded and initialized for execution. This was tested by using 
a simulated Detection Controller and a program which was executed in the other 
processor to play the role of the real-time interface function. 
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As the quarter closed, work began on a third Detection Monitor model. 
DEMON HI represents a structural rather than a functional change over pre¬ 
vious versions. The Detection Subsystem programs in DEMON III will be entirely 
core-resident at all times. In previous versions, programs were loaded from 
disk memory to perform infrequent functions such as handling I/O device errors 
and operator communication. The delays introduced by such an arrangement 
cannot be tolerated in a real-time system. Therefore an appropriate subset 
of the Disk Operating System transient functions will be made a part of the 
permanently resident nucleus. 


2.3.2.2 Event Subsystem 

Significant progress was made toward providing operating sytems tailored 
specifically to SAAC needs. These systems take the form of modifications to 
the Disk Operating System. The manner in which the modifications are provided 
allows the generation of a standard DOS, an SAAC tailored model, or anything 
in between, by specifying parametric values associated with one set of macro 
source statements. That is, only one system source copy must be maintained. 

The first step in providing suitable systems was to include the support of 
devices in the SAAC configuration which are not supported in DOS. This fulfills 
our needs for program debugging at SAAC and serves as the basis upon which 
other monitors may be developed. 

These changes to the supervisor were implemented and tested during this 
quarter, and the package was provided for daily general-purpose operation in 
place of the standard system. It includes input/output macros and error cheeking 
for the following' devices: 

1G27 Plotter 
Channel-to-Channel 
2250 Model I Display 

Experimental Display/Strip Chart Recorder 

Also, an effort was made to provide a higher level of 2250 Display support 
by moving a capability from the IBM/3G0 Operating System into the SAAC system. 

A fundamental requirement for the monitor needed to support development 
of event processing programs and to support the ultimate automatic operation is 
the capability to schedule and execute tasks automatically, as previously out¬ 
lined. The capability to dynamically schedule tasks into the background job 
queue has been designed and coded as an extensive set of modifications to the 
Disk Operating System. This scheduling program must also coordinate the 
spooling of outputs generated by two programs executing simultaneously and 
intended for the same destination (e.g., printing). These functions were under¬ 
going test as the quarter ended. An extension of this will allow a list of jobs 
to be created and managed dynamically for the foreground partition as well as 
the background partition. This work is in the design phase. 
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Other support includes the specification for a possible remote terminal 
operation which can be built into the monitor discussed above. 


2.3.3 Plotting Routines 

The SAAC machine configuration contains an on-line IBM 1627 plotter. 

The step size of this particular digital plotter is 0.01 inch. The pen can be 
moved in three directions: horizontally, vertically and diagonally. Moves 
which are not in one of these basic directions are accomplished by a combina¬ 
tion of the basic movements in such a way that the approximation is always 
within one-half the step size. 

To facilitate the use of the plotter, several subroutines were written for 
use with FORTRAN programs. Instructions for the plotter are written on mag¬ 
netic tape by these routines and later sent to the plotter by a utility program 
which may be run simultaneously with other programs. 

Three basic routines are available for line generation and annotation of 
plots. Their FORTRAN names are PLOT, SYMBOL , and NUMBER. 

The PLOT subroutine calculates the necessary plotter instructions to move 
from one coordinate to another with the pen raised or lowered. The coordinates 
are given in inches of deflection from a present origin. The algorithm to gener¬ 
ate lines which are not in one of the basic directions is an integral part of the 
routine as is the handling of all output instructions. 

The SYMBOL subroutine accomplishes the translation of alphameric infor¬ 
mation into sets of pen movements. Titles and other annotation, as well as 
symbols drawn at specific data points, can be drawn from lists of letters on 
special codes. 

Numeric data is translated to its alphameric equivalent, with the appropriate 
number of decimal places, by the NUMBER subroutine and then sent to the 
SYMBOL routine. 


2.3.4 Display Program 

Design and implementation of an initial program capability for both the 
Experimental Display and the IBM 2250 Visual Display unit were accomplished 
during this quarter. 

The 2250 display program accepts as input the seismometer level edited 
tapes and displays up to eight channels on the display. The channels are 
selectable through the 2250 f s console typewriter. Scaling is controlled through 
the functional keyboard of the 2250. Appropriate time and channel identification 
information is displayed in addition to the seismometer waveforms. 
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The Experiment Display program accepts as input the rectified and inte¬ 
grated beam tapes produced by the Detection Controller and displays the traces 
as a 32x32 matrix on the cathode ray tube experimental display. The beams 
can be displayed in a forward or reverse direction, in real time or at a fast scan¬ 
ning rate, in a single step mode in which the display is updated whenever an 
interrupt button is pushed, or in a stationary time state. The scaling and pedes¬ 
tal values which dynamically control the presentation are continuously adjustable 
from the Experimental Display Unit. 


2.3.5 Experimental Programs 

The implementation of two experimental programs—Time Delay Correlation 
and Beam Pattern Correlation—began during this quarter, and an initial capability 
was developed for both areas. Continued modification of these programs is antici¬ 
pated as experimental results evolve. They are discussed briefly in the following 
paragraphs. 

The Time Delay Correlation program was written to provide an automatic 
means for determining subarray beam steering delays. Using an iterative pro¬ 
cess, the program performs a cross-correlation between an array beam and each 
of the subarray beams or sensors in turn, adjusting each subarray beam delay 
based upon the incidence of the maximum cross-correlation. A more detailed 
description will be found in Appendix II-1 of this report. 

A Beam Pattern Correlation program was developed to scan the Detection 
Processor Short Time Average (STA) output for seismic events. The program 
accepts a reference pattern as input along with the STA data and generates a 
display tape as output as well as the printed output. The correlation is performed 
over each STA time sample of 1024 Beams to generate a set of correlation indexes. 
These indexes represent the correlations for each position of the reference pat¬ 
tern in the beam field. These values are then adjusted and placed on the display 
tape in such a way that the positions with the largest correlations will be shown as 
high intensity areas on the display. Each input STA time sample generates one 
display time sample. This work is currently in progress. 
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Section 3 


PLANS 


The following activities will be initiated during the next quarter. Discussions 
correspond to the functional disciplines. Each task is presented in sufficient 
perspective to identify its interrelation and contributions to the project. 


3.1 OPERATIONS 

Activity for the next reporting period centers upon completing plans for SAAC 
machine expansion to three computers, as indicated in Appendix 1-1. 

Plans call for machine utilization to increase with a partial second shift opera¬ 
tion scheduled throughout the period. In addition to program verification and debug 
efforts continuing at an accelerated level, two other activities will utilize major 
blocks of computer time. These activities include the generation of a modest library 
of edited LASA event tapes and an expanded effort in data analysis of these events. 

The availability of applications programs will allow investigations of the means 
of utilization, limitations, and capabilities of both the IBM 2250 CRT Display and 
the Experimental Beam Display. The evaluation of these displays is expected to 
have a major effect on the definition of requirements for system maintenance and 
operations consoles. 


3.2 SYSTEMS 

The time delay correlation program reported in Appendix II-1 appears to war¬ 
rant extension. Emphasis will be placed on experimentation for its intended pur¬ 
pose and upon applying the technique to event location. In the area of computa¬ 
tional signal processing techniques, null steering development and spatial beam 
detection will be investigated as potential system functional components. 

With respect to the time delay library process structure, the chief tasks 
ahead are the following: 

a. Specify, in detail, a maintenance procedure and schedule for checking 
the library of phase delays 
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b. Specify the conditions under which the correlation program will be 
used for updating, correcting, and enlarging the library 

c. Establish criteria for identifying events to be used in library 
updating. 

Programming of the sector interpolation formulas will be undertaken. It is 
hoped that the effectiveness of these formulas in obtaining accurate steering 
delays will be tested by applying them to several events located outside the 35 
regions. 

The magnitude-kinetic energy "density 11 relationship will be used to find 
event magnitudes for a number of events using unfiltered and filtered LASA 
beams and various window sizes. Comparisons will be made with event magni¬ 
tudes calculated by the classical method. The differences in magnitude esti¬ 
mates of each LASA subarray to that obtained on an array basis should be 
explored for different values of magnitude. 

An attempt will be made to find polynomials which fit body wave travel time 
and horizontal inverse phase velocity as functions of both epicentral event range 
and event depth. 

In addition, it appears desirable to investigate certain system instrumentation 
configurations to establish possible interim on-line operating objectives and 
requirements. 


3.3 PROGRAMMING 

The planned effort includes activity in the Detection and Event Processor 
areas, as well as in the continued development of experimental and data analysis 
programs. 

Primary emphasis in the Detection Processor area will be on extending the 
existing capability of the real-time version of the Detection Controller to operate 
with varying array configurations. This extension will provide the mechanism 
to evaluate storage and timing considerations as well as detection capability for 
different magnitude events. In the Detection Monitor area, design and imple¬ 
mentation will proceed on DEMON III, the core resident version. A version of 
the Detection Controller will be merged with the first Detection Monitor Model 
and tested by simulation. 

In the Event Processor area, the present supervisor will be extended to 
allow the experimental display to function more efficiently with the processor. 
Currently, when the experimental display program is executing, the Event 
Processor must be dedicated to this function, thereby preventing the normal 
multiprogramming mode of operation. 
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In the experimental and data analysis programming areas, efforts will be 
continued in developing and extending the experimental and IBM 2250 display 
capability, time delay correlation, noise correlation, and beam pattern correla¬ 
tion areas. 
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Appendix I 


SAAC COMPUTER EQUIPMENT 


Figure 1 depicts the configuration of SAAC machines at the end ot this 
reporting period. Changes made to the configuration sinee the last report fall 
into two categories: 

a. Additions made to implement the planned growth of the laboratory. 
Included in this category are the addition of the IBM 2250-1 Display 
Unit, the IBM 2740 Remote Printer Keyboard, the IBM 029 Inter¬ 
preting Keypunch, and the expanded capability Experimental 
Display. 

b. Additions and/or reconfigurations undertaken to better adapt the 
system to its evolving utilization. Because of the heavy emphasis 
on program evaluation and debug, independent operation of each 
System/360 Model 40 became increasingly desirable. As a result 
a 2841-1 Disk Control Unit and 2311-1 Disk Unit were added to 

the Detection Processor so that each processor could operate under 
an independent Disk Operating System. 

Integration of the Experimental Display with the Event Processor has 
released a 2403-2 Tape Drive and Control Unit. These were connected to the 
Detection Processor Channel No. 2 to provide a tape dump for print-out. The 
resulting tape is printed by the Event Processor on the 1403-2 line printer 
operating in a foreground mode and overlapped with normal Event Processor 
operation. 

Figure 2 shows the machine system developments planned during the next 
reporting period. The addition of the 1403-2 Line Printer to the 40G will 
eliminate the need for the 2403 Tape Control and Tape Unit and allow it to be 
deleted from the system. The replacement of the 2520 by the 2540 Reader-Punch 
will improve the efficiency of card operations. 

Figure 3 illustrates the planned Auxiliary System configuration. 
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Figu e 1 Co rent Eqoipmen* Configuration 
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Figu e 2 Planned Equipment Configuration 
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Figure 3. Auxiliary Equipment Configuration 





























Appendix II 


SIGNAL PROCESSING 


In this appendix, the technique of automatic time delay generation and the esti¬ 
mated event location accuracy of large arrays are presented. The results of a 
study on the use of fidelity optimum processing against one of two simultaneously 
occurring events appear promising and warrant both study continuation and initial 
experimentation. The accuracy of fixed point arithmetic execution and the esti¬ 
mated execution speed improvement obtained by microprogramming the Fast 
Fourier Transform is discussed. The results of studies in array optimization are 
also included. 


II. 1 TIME DELAY CORRELATION 


II. 1.1 Introduction 

The time delay correlation procedure described herein was suggested as a 
topic of investigation because it was felt that the large volume of time delay data 
anticipated during system operations could not be efficiently and precisely handled 
manually. The results of this study and especially its implications on event pro¬ 
cessing warrant further priority considerations. In the present section, the method 
and some experimental results are given and the effects of incorporating wavefront 
analysis are also discussed. The potential effects of the findings on system per¬ 
formance and recommendations are given in the conclusions. 


II.1.2 Procedure 

Depicted in Figure 4 is the formulation of the time delay correlation procedure. 
It employs plane wave delays with regional corrections (see Appendix III-2) as 

initial values for the time delays (Tj_, t 2 ,-, r n ). The Laplace notation is used to 

signify the delay process from which the beam [B(t)] is formed. Each input pc 1 (t), 

X 2 (t),-,X (t)] is correlated with the beam. The time (f n ) at which the peak value 

of the correlation function [R n (r)] occurs is selected and used to estimate the cor¬ 
rection to the initial time delay values. The quantity a has been introduced to 
provide a process stability control parameter (i.e. when a 1, full value correc¬ 
tions are applied). In the process programming, operational coding notation and 
floating point arithmetic have been employed. 
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Figure 4. Time Delay Correlation 
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Once initialized, the process is allowed to iterate K times. As soon as 
f n (K^l) is not significantly different from r n (K), it is conjectured that the beam 
power is near maximum and the time delays are optimal. 

For the limited testing performed to date, the procedure has been found to be 
susceptible to false nulls. As expected, it has been observed that the correlation 
function is sinusoidal and that suboptimal convergence occasionally occurs at an 
integral number of periods away from maximum. The system’s narrow band 
response can be employed to remove this ambiguity by using wave front analysis 
techniques to test the observed nulls for credibility. Since the false nulls are large 
with respect to wave front corrugations, the procedure is expected to yield satis¬ 
factory results. 


II.1.3 Experimental Results 

To demonstrate technique, a limited sample of certain experimental results 
is given. Although use of subarray beams would no doubt yield better results, 
single seismometers were employed for computational convenience. Depicted in 
the upper portion of Figure 5 is a well-behaved element trace (a) and its beam 
correlation function (b). A poor element signal (c) and its beam correlation func¬ 
tion (d) are shown by the lower pair. Note that the sensitivity of (b) is much 
greater than (d) and, hence, for this data, the confidence of estimating the time 
delay from (b) is significantly greater than that for (d). 

The effect of signal-to-noise ratio was also investigated by applying the cor¬ 
relation procedure to scaled Longshot data for which the average input signal-to- 
noise ratio was significantly reduced. The average of the magnitude of the differences 
between the correlation time delays and the time delays determined from the 
unsealed Longshot traces was computed. Denoted by A in Figure 6 are the results 
for three scaled sets of data with input signal-to-noise ratios of lOdB, 5dR and 
-1.5dB. The effects of repeating correlation after employing wavefront corrections 
are indicated by □ , in Figure 6. 

Typical waveforms of the scaled and unsealed traces are shown in Figure 7. 

For the latter, the signal-to-noise ratio is almost 35dB and the time delay can be 
determined with little difficulty. However, when the signal-to-noise ratio is reduced 
by 36dB so that the input signal-to-noise ratio is negative, and even when the 
traces are appropriately time-phased as is the case in Figure 7, it is obvious 
that a manual time delay measurement is impractical. 


II.1.4 Conclusions 


The importance of accurate time delays cannot be over emphasized. First, a 
loss in beam gain is experienced^ anc ] the resulting beam defocusing will degrade 
array/processing performance. Second, event location errors (see Appendix II. 2) 
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Figure 6 . Correlation Estimation of Phase Delays 
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Figure 7 Seismorre'er Record cf Longshot 

















are strongly dependent on the time delays. Third, the maintenance of both the time 
delay and location conversion libraries in a form suitable for updating implies an 
automatic processing procedure for both time delay determination and use. 

The method discussed above provides the desired results. The sparse experi¬ 
mental results available to date warrant consideration of the technique as an integral 
part of event processing. 


II.2 EVENT LOCATION ACCURACY 


II.2.1 Introduction 


The error generated when estimating the location of a seismic event received 
on a large array can be attributed, in part, to the presence of seismic noise at the 
time of arrival of the event signal (i.e., to the fact that the signal-to-noise ratio 
is finite), and, in part, to imperfect "calibration." In the present context the term 
"calibration" refers to the method and data used to convert wavefront measurements 


(such as the reciprocal wave velocities u 


, u ) to a geographic location for the event. 

x y 


In Section 2, an idealized model is used to examine the effect of seismic noise 
on the accuracy of the estimated event location. In Section 3, an attempt is made to 
incorporate the errors resulting from imperfect calibration. 


II.2.2 The Effect of Seismic Noise on Estimated Event Location 

Because of the presence of seismic noise at the time of arrival of a seismic 
event at the array site, errors are generated in the estimated values of the recipro¬ 
cal phase speeds (u x , Uy) of the event wavefront. Hence, even with perfect "calibra¬ 
tion", errors in the geographic event location will result. These errors tend to 
increase with a decreasing signal-to-noise ratio. 

To assess the effect of seismic noise in a quantitative manner, we consider an 
idealized large array model, wherein array beams are formed from subarfay beams 
for which the following behavior is postulated: 

a. Seismic noise received at any one subarray is statistically independent from 
that received at any other sub array 

b. The signal waveforms received at all subarrays are identical except for 
(frequency independent) scale factors and for travel time differences 

c. The wavefront shape is accurately known and the time delay for the k-th 
subarray can be represented as t^ + x k u x f Yk u y» where t\ K is the (accurately 
known) station correction for the k-th subarray and (x^, y^) are the posi¬ 
tion coordinates of the center of that subarray. 
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The first assumption is in reasonable accord with present experience. The 
second assumption is not always fulfilled: substantial decorrelation has been 
observed for signal waveforms received on different subarrays. However, by and 
large, the correlation is satisfactory for a few seconds following the signal onset. 
The third assumption postulates perfect "calibration" of the wavefront shape. 

This assumption is realistic if the seismic event originated in close geographic 
proximity of previously recorded strong seismic events. Where such previously 
received events are lacking, the assumption is not always valid. 

With the postulates, the estimated reciprocal phase speeds (u x , u v ) approximate 
Gaussian variables with the true quantities (u x , Uy) as mean values. Furthermore, 
assuming application of (frequency independent) weighting of the subarray beams to 
achieve maximum array beam output signal-to-noise ratio, the contours of constant 
probability density for (u x , u, r ) coincide with the constant-loss contours of the beam 
pattern centered at (u x , Uy). ^Such weighting is desirable, for example, when the sub¬ 
arrays do not all have the 7 same number of seismometers. 

Denoting by p the probability for (u x , Uy) to fall outside the beam contour cor¬ 
responding to a loss of L dB for a center frequency of f Hz, the following approxima¬ 
tion holds: 



(1) 


where ctq is the standard deviation of the timing error for the signal waveform in 
the LASA beam output and is given by: 



( 2 ) 


Here S/N is the output signal-to-noise ratio of the array beam steered to the event 
location (u x , Uy), while WT is the time-bandwidth product (for LASA, WT ~ 1 ). 

The above formulae provide minimum error bounds and tend to under estimate the 
errors in (u x , Uy) particularly for low S/N. 

The preceding results conform to the intuitively obvious fact that the location 
error (u x , Gy) can be minimized by positioning the subarrays in such a manner 
that the resulting beam pattern has minimal main-lobe width. For the particular 
case where the subarrays are constrained to be located in a circular area of 
radius R, they should be positioned on the circumference. In that case formula 
(1) reduces to: 


p^exp -( 7 rfrR) 2 • ^ • 2WT 


(3) 


where r is the distance between the estimated location (u x , Uy) and the true location 
(u x , Uy). The root-mean-square of r is given by: 



( 4 ) 
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The error, r, is measured in see/Km and is related to the (u x , Uy) plane. To 
translate it as a geographic location error, it must be multiplied by a sensitivity 
factor C which generally depends on the range from the event location to array as 
well as on the azimuth from the true event location to estimated event location. 

After averaging over this azimuth, the factor C varies from about 1.2 x 10^ to 2.0 x 
105 Km^/see. Thus, C = 1.6 x 10^ Km^/see appears to be a satisfactory compro¬ 
mise particularly suitable for ranges from the event location to the array in exeess 
of 65°. 

As an example, consider the case where f = 1 Hz, WT = 1. Then, the root- 
mean-square geographic error is given by: 


„„„ , . 0.36 x 10 N ' T , 

RMS geographic error « -—- J — Km, (5) 

with the array radius R being measured in Km. This result is depicted graphically 
in Figure 8. 

We emphasize that minimization of the main-lobe width may not be a praetieal 
design criterion since it can lead to a poor side-lobe structure. 

The expressions (3) and (4) presuppose that all subarrays show the same output 
signal-to-noise ratio. When this is not the case, the value of S/N in (3) and (4) must 
be derated. For example, when all subarrays have the same number of seismometers, 
except for one subarray whieh has a larger number, the value of S/N must be derated 
by less than 0.4 dB, provided that the large subarray has no more than half of all the 
seismometers in the array and provided that the subarrays are suitably positioned 
relative to eaeh other. 


II.2.3 The Effeet of Imperfect n Calibration M on Location Error 

The considerations in Section 2 have assumed perfect ’Calibration" (i.e., have 
disregarded the errors arising when converting the estimated reciprocal phase 
speeds (u x , fty) to geographic event location). These errors exist even when the 
event depth is known or postulated. 

The main eause for "calibration" errors lies in the fact that strong seismic 
events, needed to establish the data base for the ’Calibration" technique, are rather 
rare; favor specific geographic areas to the virtual exclusion of other large regions; 
have finite signal-to-noise ratio; and are sometimes difficult to pinpoint in depth 
and geographic location. These eircumstanees not only impede precise determina¬ 
tion of the relative station corrections required for beamforming, but also hamper 
the establishment of an accurate conversion method from (u r , u ) to geographic 
location. X V 


Although extensive data analyses have been undertaken to obtain the relative 
station corrections for LASA, no explicit study has yet been made to calibrate a 


conversion method from (u 

' V 


u ) to geographic location. 


Presently, ray tracing 
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techniques, based on the Jeffreys-Bullen tables and others, are used for this 
conversion, but no systematic effort has been undertaken to compare the calculated 
event locations with those obtained from other observations (such as from a net¬ 
work of seismic stations) and to establish corrections for possible bias in the event 
locations calculated from LASA data. In particular, little information appears to 
exist regarding the event location accuracy presently achieved by LASA. 

One may conjecture that the n ealibration M error will eonsist of two portions: 
one part being proportional to the reciprocal of the array diameter, 2R, and 
refleeting the effect of the finite array aperture; the second part being independent 
of aperture size and reflecting the location inaccuracies for seismic events used in 
the data base. We shall assume that both portions ean be looked upon as random 
errors when considered on a worldwide basis. Hence, combining these errors in 
RMS fashion with the error presented by formula (4), we obtain the following expres¬ 
sion for the total geographic RMS location error: 



A 


+ B Km, 


( 6 ) 


+ 


R 


where the constants A and B will depend on the '’calibration” accuracy. 

Although in the initial phases of the operation of an array, large "calibration" 
errors may exist (reflected in large values of A and B in formula (6)), as time pro¬ 
gresses and the database broadens, the "calibration" errors should decrease and 
the values of A and B should become smaller. Also, A and Bwill tend to be less 
for geographic areas where seismic events are relatively prevalent than for areas 
where they are rare. 

The results of formula (6) have been plotted in Figure 9 using the same data 
as in Section 2 (i.e., taking f=l Hz, WT=1, and Ol.GxlO^ Km^/sec) and choosing 
A=10* Km^ and B=0. This choice for A and B implies a calibration error of 100 
Km for an array radius of 100 Km and has been somewhat arbitrary. 

It is interesting to compare this error with the location error which would be 
incurred when the sole contributing eause is assumed to be the fact that the sampling 
rate is finite, thus limiting the accuracies with which the time delays for array 
beamforming ean be implemented. This location error ean be calculated from 
formula (1) provided that the standard deviation be defined by: 

2 _ (AT) 2 

a 0 12. K 

(based on the assumption that timing errors due to time sampling are uniformly 
distributed over the sampling interval AT and are uneorrelated between subarrays). 
Here, K is the number of subarrays. For an array whose subarrays are of equal 
size and are positioned on a circle with radius R, the RMS error is found to be 
2 (7 0 C/R. For a sampling rate of ten samples per second (AT 0.1 sec), K 21 and 
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Figure 9. Location Error with Imperfect "Calibration" 
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5 2 

C 2x10 Km see, one finds a location error of 20 Km when It 100 Km. Thus, we 
assumed for Figure 9 a "calibration" error five times larger than the error caused 
by the sampling rate. For LASA, formula 1 leads to an It MS error of 40 Km, 
caused by sampling rate. 


II.3 FIDELITY OPTIMUM PROCESSING 


II.3.1 Introduction and Summary 

This appendix contains the results of a theoretical study to estimate the ability 
of fidelity optimum processing to discriminate against one of two simultaneous 
events. It compares the signal-to-noise ratio of conventional processing and fidelity 
optimum processing when the noise consists, in part, of a second, interfering seis¬ 
mic event. Antenna patterns of the resulting optimum processor are also obtained. 

Specifically, the noise is assumed to be the sum of an uncorrelated background 
and an interfering event. The background is assumed to have the same power spec¬ 
trum at each seismometer. The interfering event is also assumed to be a second 
order stochastic process with a power spectrum identical to that of the background 
noise. 

In this situation, the gain from optimum processing is a function of the array 
configuration, the radius of the array, the ratio of the power of the interfering 
event to the background noise power, and the magnitude of the difference in inverse 
phase velocity between the event of interest and the interfering event. Three array 
configurations were considered: The LASA configuration, 21 seismometers equally 
spaced on a circle and five seismometers equally spaced on a circle. If we let R 
be the radius of the array and |Au | be the magnitude of the difference in inverse phase 
velocity between the two events, then the gain is a continuous function of R | Au| and 
we determined this function over the range of interest. Ratios of interfering event 
noise to background noise varying from 0.05 to 100 were considered. 

Substantial gains were found to occur for r|au | > 0.2 for the LASA configuration 
and for R| Au| >0.15 for the other configurations. In the above ranges of R |au| and 
for large ratios of interfering event noise power to background noise power, the 
optimum processor essentially eliminated the interfering event. 


II.3.2 The Power Output of a Linear Processor 

We suppose we have K seismometers and designate their sampled outputs by 
X.(t),..., X^(t). We assume: 

X k (t) = S(t) + N k (t), k = 1.K; (8) 

where each seismometer output is the sum of a signal, which we are interested in 
estimating, and noise. The signals are assumed to be lined up or in phase from 
seismometer to seismometer. We suppose the noises, N^(t), to form a multi- 
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variate second order process with a spectral matrix P^,(f). That is, 


E{N k (t) N k ,(t+r)} = <^ k ,(r) 


and 

* (9) 

P kk' (f) At £, \k' (T) e ^ ’ 

7*=_oo 

where At = sampling interval. 


Now we will consider the case where a linear operator with impulse response 
0 k (s) and transfer function A k (f) is applied to the kth seismometer output for k=l, 

.. .,K. These filter outputs are summed to obtain the output of the array, Y(t). This 
is illustrated in Figure 10. 


0 k (s) and A k (f) are related by: 


= ^ A t (t)e 2l)fs df, 


V s > 


-w 


and 


A fc (f) = At 0 k (s) e 


- 2 7Tj 1S 


S =-00 


where w - 


2 At 


( 10 ) 


We will be interested in the output noise power of an array with such processing. 
The output noise power spectrum (see Figure 11) is given by: 

1 V) A k ,(f) p kk ,(f). (11) 

k^; 

By A k (f), we mean the complex conjugate of A k (f). 

Hence, the output noise power is given by: 
w 

I I A k' f > A k'« t > p kk' (f > df - < 12 > 

" w k,k' 

II.3.3 Form of the Fidelity Optimum Processor 

In fidelity optimum processing, the goal is to minimize the output noise power 

subject to the constraint that the signal is undistorted. Now, the transfer function 

K 

seen by the signal is Z A, (f); hence, the goal of fidelity optimum processing is to 

k=l K 

minimize: 
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Figure 10. Linear Processor 


30 











0 

5 

10 

15 

20 

25 

30 



4 5 

Frequency (Hz) 


Figure 11. Noise Spectrum Used 


31 
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S 1 V 


c o A k- (f > p kk . <'> 

-w k,k' 

subject to the condition: 


(i:t) 




(f) = i. 


(14) 


7 2 8 

It can be shown (see Capon or Vanderkulk ’ ) that the Aj_(f) satisfying the above 
are given by: 


I p kk’ 


<f) 


A k (f) - 


£ P ki> 

k,k' 


(15) 


where P. , , (f) is the k,k' element of the inverse of the KxK matrix P, , . (f). 
k,k' ' kk'' ' 

II.3.4 Inverting the Spectral Matrix 

Pkki (f) for fixed f is a matrix of complex numbers. If we let: 

P kk- < f > - C kk' ([ > + J «V <*> (16) 

and 

p kk' (f > = D kk'< f > + J E kk'< f >- 

then it is easily seen that the following relationship holds between real 2K x 2K matrices 




-Q 


D - E \ 




Q 




E 


D 


(17) 


/ 


,-l 


This relationship can be used.to find P, , , (f) with a real matrix inversion routine. 
We can also use the above form to show that 

P kk' (f) = p Do. (*) and ’ hence ’ that : 


k'k 
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.-1 


kk' | kk' ( l_, kk' 

k,k' k,k' ^ J k,k' 

C -Q \ 

Consider first the matrix | I . Since C is symmetric and Q is skew 

Q Cj 


( 18 ) 


D 


symmetric, this matrix is symmetric. Hence, its inverse 


, is sym- 


D 


metric which means that E = -E. 

By E 7 we mean the transpose of E. This is equivalent to P (f) = P^q. (f) • 


II.3.5 A Convenient Form for the Output Noise Power with Optimum Processing 
We saw earlier that the output noise power in optimum processing was given 


by: 


w 


^ ‘ ^ A k( f ) A k' (f) P kk' (f) df ' 


(19) 


-w kk' 

This is an integral over a quadratic form which can be written as: 

I V> A k' < f > P kk- < f > = A T P A. 
kk' 



( 20 ) 


and P is the cross spectral matrix. We have suppressed f for simplicity. However, 
we also saw that: 


> p ki> 

T' 


AJf) = 


p kk- ® 


= P _1 A , 


( 21 ) 


kk' 


where 


A = 


P kk' ( f ) 


(\\ 

i / 


kk' 
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T -1 -1 

/"«w» * X /-«W» T _ I rvj T TT 

Hence, A PA=AP P P A = A P A 

-1 ^~ 1 ~^ ~ ~T ^ -1 

However, since P^, (f) = P ttt (f), we have A = A and P = P 


kk' 


Hence, finally : 


~ T T -1 kk 1 

A P A = A P A = — 


I P kk' (f| 


P kk' (t > 


kk' 


(T p kk '< f >'' _1 

\ kk 1 


( 22 ) 


or 


A P A = 


This yields • 


-1 \ 
Re { p kk - ( f » 


kk' 


-1 


/ 


W 

output noise power = II A k ( f ) A k .( f ) p kk .( f ) 


-w kk . 


(23) 


(24) 


or 


w 


-II Re { p kk-i (t >i 


-w 


kk' 


df. 


(25) 


II.3.6 The Output Noise Power in Conventional Processing 

In conventional processing the output noise power is given by: 

U/ 

(1/K 2 ) ^ Re{P kk ,(f)} df 


I 


-w 


kk' 


II.3.7 The Special Case of Noise Consisting of a Single Seismic Disturbance Plus 
an Uncorrelated Background 

Suppose the noise n k (t) to be the sum of two noises n^ k (t) and n k (t). Further 
suppose n c k (t) to be a seismic disturbance with delays t . ,r relative to the 
delays of S(t). Then: K 


n c f kW = 


R c^- T k 


). 


( 20 ) 


34 









Further suppose that n c (t) is a second order process with power spectrum (f) 
and autocovariance function (p ( T )- Then: 

E{ n ck (t) n ck ( t + T )l = E ( n c ( t_T k ) n c ( t_T kt + T )^ (27) 


= <P(t+ T k - r k ,) 

Taking the Fourier transform we obtain P c (f) e J ' k k !/ as the contribution 
to the cross spectrum P^ k , (f). 

Now njj k (t) is assumed to be a background noise which is uncorrelated from 
seismometer to seismometer and uncorrelated with the seismic disturbance. If 
it has a power spectrum P^ (f) then the cross spectra are given by 

P kk . (t) = P b ( f ) <5 ( k > k ’) + p c ( f ) e 2 ?rjf(T k" T k’ ) ; (28) 

or if we let P (f) / P (f) = M(f) } 

0 u 

P kk' (f) = p b ( f ) { 6 ( k - k ’) + M(f) e 2?rjf(T k " V}, (29) 

8 2 ‘j 

where <5 (kk f ) is the Kronecker delta. From Vanderkulk and Kelly and Levin * we have 



(f) 


p b ( f ) 


{6 (kk’) - M (f) e 27r j f(T k " T k !) / (1 + KM (f))} . 


( 30 ) 


Using the above we can write down expressions for the transfer functions A^(f) 
and the noise power output for optimal and conventional processing. The 
transfer function is: 


A k ( f ) 


Ml 

l+KM (f) 

M (f) 

l+KM (f). 



e - 2?r j fT k' 



2 7Tj f (\ 


T w) 


kk' 


(31) 


The output noise 



(optimum processing) is: 


M(l) 

l+KM(f) 


I 


e 2 7T j f ( T k - V 



df, 


(32) 
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and the output noise power (conventional processing) is: 


w 


f P b (f) {K+M (f) ^ e 27r ^ f(T k " V} / K 2 df. 


( 33 ) 


kk' 


The calculations can be simplified by noticing that: 


2 


e 2 ') f < T k "V 



27t jf x 1 2 

e kl 


(34) 


kk' 


k 


II.3.8 Determination of Delays 

We suppose K seismometers are located by the vectors p^ with polar coordinates 
(Rr k , T^). R is a parameter which controls the overall radius of the array. We let 
u”l be the inverse phase velocity^ofJ,he event of interest,"^ be the inverse phase veloci¬ 
ty of the interfering event and Au=u 1 -U 2 * Further we let Au have polar coordinates 
( |au |, <p). With these definitions the delays, t, , of the interfering event with respect 
to the event of interest are given by 


T k = Au . p k 


(35) 


= R |Au | r k cos (y k - <p) 

In our investigation the will be normalized so that their maximum is approximately 
equal to one. Hence, the delays are a function of R |Au( and <p. 

II.3.9 The Antenna Pattern 

For the class of linearly processed arrays under consideration, the antenna 
gain in the direction with delays r^, .. ., is : 


I A k < f > 



2 


antenna gain (f, t^, . . . , r.) 


(36) 


k 
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This is a power gain, and if we have an input from this direction with power 
spectrum P(f), the power output will be: 


w 


( 37 ) 


I p,f > I 



-w 


k 


II.3.10 Specific Calculations 

We now give the results of some calculations to determine the ability of optimum 
processing to discriminate against the noise described in Section II.3.7 (i.e., noise con¬ 
sisting of a single seismic disturbance plus an uncorrelated background). 

In these calculation, we assumed P c (f) and P|j(f) to have the same shape: i.e., 

P c (f) = P(f)=MPb(f)- We took as the shape of P(f) a smooth approximation to the 
spectrum of the output of an individual seismometer for the Kamchatka earthquake 
(April 8, 1966) as given in Figure 3-43, Page 3-130 of Ref. 4. 

As the ratio M of the spectra P c (f) and P| 3 (f),we took values from 0.05 to 100. We 
considered three configurations of seismometers. The first was that for the 21 
subarrays of the present I.ASA. The r^ andy^ are given in Table 1. The second 
was an array of 21 seismometers located equally spaced on the circumference of 
a circle. The third was an array of five seismometers located equally spaced on the 
circumference of a circle. 

For these configurations and for a wide range of directions of arrival of the 
unwanted seismic event, we obtained the noise power output for optimum processing 
and for conventional processing. We also computed the ratio of the conventional 
noise power output to the optimum noise power output. In addition, we calculated 
antenna patterns for both a sample case of optimum processing and conventional 
processing. 

In these calculations, we assumed that the array was pointed in the direction 
corresponding to all delays equal to zero. The noise event will then have a set of 
delays, not all equal to zero, which we saw in Section II.3.8 are functions of It |Au| and 
0. where R is the radius of the array and (Au, $) are polar coordinates of the difference 
in inverse phase velocity of the event of interest and the interfering event. We did 
the calculations for various values of <t> and generally for a range of It |Au| >_ 0. The 
results are plotted against R |au | with <p fixed. 

II.3.10.1 Results of the I.ASA Configuration 

In Figure 12, we have plotted the noise power output for conventional processing 
and, in Figure 13, the noise power output for optimum processing as a function of 
the direction of the interfering seismic event. This is done for values of M from 
0.05 to 100. Specifically, the noise spectral matrix is assumed to be: 
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Table 1 


LAS A COORDINATES 


Instrument Group 

Normalized Position 

r k 


AO 

0.000 


B1 

0.123 

55.2° 

B2 

0.007 

141.7° 

B3 

0.082 

246.1° 

B4 

0.090 

346.7° 

Cl 

0.183 

23.8° 

C2 

0.163 

97.8° 

C3 

0.131 

192.5° 

C4 

0.128 

293.7° 

D1 

0.305 

56.8° 

D2 

0.263 

141.9° 

D3 

0.252 

231.9° 

D4 

0.308 

336.1° 

El 

0.542 

13.6° 

E2 

0.686 

106.9° 

E3 

0.607 

183.3° 

E4 

0.538 

278.3° 

FI 

1.093 

46.4° 

F2 

1.037 

147.4° 

F3 

1.036 

219.4° 

F4 

0.973 

235.4° 
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Figure 12. Output Noise Power, Conventional Processing, 
LASA 
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Figure 13. Output Noise Power, Optimum Processing, LASA 
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P kk ,(f) = P(f) {6 (kk') + Me 


2ri f (r k - r k ,) } 


(38) 


where the r k and T k , are determined by the direction of the interfering seismic 
event relative to the event of interest as described in Section II.3.8. We have assumed 
0 = 0 (i.e., all events are arriving from the same relative azimuthal direction) and 
varied the r|Au| for different values of M. 

In Figure 14, we have plotted the dB gain in noise rejection of optimum process¬ 
ing over conventional processing for the LASA configuration and values of M from 
1 to 100. Notice that substantial gains occur for r|au|>0.2. 

In Figure 15, we have plotted the antenna gain of the optimum processor for an 
event with R|au| = 0.3, 0=0 and M = 10. Cuts radially out from the center are 
plotted for = 0, n/8, ..., n . Notice the sharp null in the direction of the interfering 
event. Figure 16 gives several radial cuts of the conventional antenna pattern. 

II.3.10.2 Results for 21 Seismometers Located Equally spaced on a Circle 

Figure 17 contains the noise power output for conventional processing and 
Figure 18 the noise output power for optimal processing for the case for 21 seismo¬ 
meters equally spaced on a circle for M from 0.05 to 10. In Figure 19, we have 
plotted the dB gain in noise rejection from optimum processing for M from 1 to 10. 
There is a significant improvement over the 21 seismometer configuration of the 
LASA array. Figure 20 contains the radial cut of the optimum processor antenna 
pattern through the interfering event for the same interfering event depicted in 
Figure 15. 


II.3.10.3 Results for Five Seismometers Located Equally Spaced on a Circle 

Figures 21, 22, 23, and 24 give some results for a configuration of five seismo¬ 
meters located equally spaced on a circle. Figure 23 gives the antenna pattern of 
the optimum processor again for the same interfering event as in Figures 15 and 20. 
Figure 24 gives the antenna pattern of the conventional processor for 0<0<7r /5. The 
remainder of this pattern can be obtained by symmetry. 

II.3.11 Optimum Processing Consisting of a Single Weighting of Each Seismometer 


Output Before Addition 


Suppose we wish to restrict our processing to simple weighted sum of the seis¬ 
mometer outputs. That is, we wish our output to be of the form: 



(.-!<)) 
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Figure 14, dB Gain Due to Optimum Processing, LASA Configurati 
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Figure 15. Antenna Pattern 

Optimum Processing 
Interfering Event 
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Figure 16. Antenna Pattern 

Conventional Processing 
LASA Configuration 
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Figure 17. Output Noise Power 

Conventional Processing 
21 Seismometers on a Circle 
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Output Noise Power (dB) 
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Figure 18. Output Noise Power 
Optimum Processing 
21 Seismometers on a Circle 
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Figure 19. dB Gain Due to Optimum Processing 
21 Seismometers on a Circle 
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Figure 20. Antenna Pattern, Optimum Processing 
21 Seismometers on a Circle 
Interfering Event 
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Figure 21. Output Noise Power 

5 Seismometers on a Circle 
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Figure 24. Antenna Gain 

Conventional Processing 
5 Seismometers on a Circle 
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with the fidelity constraint: 


l 



(40) 


This is the same as restricting the (f) of sections II.3.2 and II.3.3 to be of the form 
Ak(f) = cfc. With this restriction the output noise power is given by: 

w 


output noise power 


I I c k V p kk' <*> df 

-w k,k' 


(41) 


w w 

^ C k C k' P kk' (f) df 


kk' 


-w 


w 


- y °k c k' y Re { p kk' <f >) df - 


kk' 


-w 


(42) 


(43) 


Now if we let: 
w 

$ Re < p kk- (t)l df = “kk- <44) 

-w 


Then we wish to minimize 


I 


kk' 


C k C k' 


II 

kk' 


(45) 


subject to the constraint: 


I 



(46) 
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This is the same abstract problem described in Section II.3.3 and the optimum set of 
c^ are given by: 




(47) 


Further as in Section It.3.5: 
output noise power 


■ I 


C k V 11 


k,k' 


kk' 




-1 


(48) 


(49) 


and from Section II.3.9 the antenna gain is: 
antenna gain (opt. processing) = 

k 


c k e- 2l J fT k| 2 . 


(50) 


For the case of present interest, a noise consisting of an interfering event 
with spectrum P (f) and delays r, and a uncorrelated background with spectrum 
(f), we have C 

v - 1 p kk'' f » 

= 6 ( k,k') J P b (f) df + j* P c (f) cos 27Tf(T k - T Rt ) df ; 
or if P b (f) = P(f) = P c (f) /M, 


(51) 


n kk , = 6 (k,k') J P(f) df + M J P(f) cos 27rf(r k - r^) 


df. 


(52) 


We made some calculations comparing the gain in noise rejection using this 
simple processing with the gain from complete optimum processing. Figures 25 
and 26 give the noise power output for the LASA configuration and for 21 seismo¬ 
meters on a circle for single multiplication optimum processing. There is con¬ 
siderable improvement over conventional processing but substantially less than 
with complete optimum processing. For five seismometers on a circle, a sin¬ 
gle calculation for the value M=5 showed single multiplication optimum process¬ 
ing to give almost no improvement over conventional processing. 
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Output Noise Power (dB) 
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Figu r e 25 


Output Noise Power, Single Multiplication 
Optimum Processing LASA Configuration 
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II.4. FIXED POINT FFT EXECUTION 


II.4.1 Introduction 


II.4.1.1 Summary 

This appendix contains the results of an effort directed at estimating the impact 
of microcoding on the performance of the fast Fourier transform algorithm (FFT) 
on the IBM System/360 Model 40. There were two major aspects to this effort; 
estimating the speed of a fixed point. Model 40, microcode and estimating the error 
in the fixed point calculation. The microcode and its speed are discussed in detail 
in Section II.4.3. 

If N is the number of complex points to be transformed (N is assumed to be a 
power of two), then the following estimate is obtained for the speed of a fixed point 
microcode operating on 16 bit numbers; 

Transform Time ~ 52 N log i;> N ji secs. 

This is approximately ten times faster than both floating and fixed point machine 
language codes for the algorithm on the Model 40. 

The accuracy of the fixed point calculation is discussed, in detail, in Section EL4.2. 
The fixed point calculation we considered, both for microcoding and in the accuracy 
estimation, involved keeping the array properly scaled by rescaling at any stage in 
which there occurred an overflow. With this scheme, the error is a function of the 
numbers and timing of these overflows. However, an upper bound to the error can 
be obtained by assuming an overflow and rescaling at each stage of the calculation. 
Again, if we assume we are transforming N = 2^ complex points, then the upper 
bound of the ratio of the root mean square error to the root mean square of the 
resulting transform increases as the square root of N or by vlT at each of the M 
stages of the calculation. Specifically for the Model 40 microcode which operates 
on 16 bit numbers, we have: 


RMS (error) 
RMS (result) - 


2 (M+3)/2 -15 |U 3 , 

RMS (initial array) 


This is well within the accuracy requirements of the FFT applications envisaged for 
LASA data. This point is discussed, in detail, in subsection EI.2.4.5. 


II.4.1.2 The Finite Fourier Transform 

If X(j) j = 0, 1, . . . , N-l is a sequence of complex numbers then the finite 
Fourier transform of X(j) is the sequence: 
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N-l 

A(m) = (1/N) V X(j)e" 27rljm/N 


m=0, 1, . . . , N-l. 


(53) 


J=° 


The inverse transform is 

X(i) = £ A,m)e 21I ^ m//N 

m=0 


(54) 


In both of the above equations i = (-1) . We will be considering a fixed point 

calculation of these transforms using the Fast Fourier Transform algorithm 9 , 10 . 
In connection with Equation (53), we will actually obtain N A (m) from X (j). We will 
then give this result with N“1 as part of an overall scale factor. Considering the 
calculation of N A (m) from X (j) or X (j) from A (m) Parseval's theorem states: 


N-l N-l 

y N)i 2 = n > 

1 A(m) | 2 

j-o m=o 

N-l 

N-l 

or \ | NA(m) | 2 = N 

I lx «> 

L — i 

m=o 

j~° 


(55) 


and we see that the mean square value of the result is N times the mean square 
value of the initial sequence. This fact will be used below. 


II.4.1.3 The Inner Loop of the Fast Fourier Transform Algorithm; Step by Step 
Scaling 

The inner loop of the FFT algorithm operates on two complex numbers from 
the array. It takes these two numbers and produces two new complex numbers 
which replace the original ones in the array. Let X m (i) and X m (j) be the original 
complex numbers. Then, the new pairX m+ ^(i), ^ m + ^(j) are given by: 


X ^ (i) = X (i) + X (j)W 
m+1' ’ m m u/ 

X ^ (j) = X (i) - X (j)\V 
m+1 ' J/ m' ’ m 
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where W is a complex root of unity. If we write these equations out in terms of 
their real and imaginary parts, we get: 

Re{X m+ i(i)} = Re{X m (i)} + Re{X m (i)} Re{W} - Im{X m (i)} Im{w} 

(57) 

Im ^ X m+l (j) ^ = Im (X m (i)} " Re ^ X m (^ IrT W " Im ^- X m (j^ Re M 

At each stage the algorithm goe^ihrough the entire array of N numbers in this 
fashion, two at a time, If N = 2 , then the number of such stages in the computation 

is M. 

As we move from stage to stage through the calculation, the magnitudes of the 
numbers in the array generally increase which means that the array can be kept 
properly scaled by overflow tests and right shifts. Consider first the root mean 
square of the complex numbers. From equations (56) we have 



Hence, in the root mean square sense, the numbers (both real and complex) are 
increasing by V~2 at each stage. Consider next the maximum modulus of the com¬ 
plex numbers. From equations (56) one can easily show that 


mux {| x m (i)|,| x m (j)|} < max {|X m+1 (i)|, |X m+1 (j)|} <. 2 max {| X m (i) | , |x m <j)|} 

(59) 


Hence, the maximum modulus of the array of complex numbers is non-decreasing. 

Both of the above results suggest that the fixed point calculation be done as 
follows; the entire original array and the W f s are scaled with the numbers moved as 
far to the left as possible; after the multiplication, the most significant bits are 
retained; after the addition, there is a test for overflow; if overflow is detected, the 
entire array is shifted right one bit and the calculation is continued. Equation (59) 
shows that the maximum modulus of the complex numbers cannot increase by more 
than a factor of 2. Hence, if scaling were controlled by this modulus, there would 
be no more than one shift per stage. However, the method we have outlined above 
would control the scaling on the basis of overflow of the calculation of the real and 
imaginary parts. It can be shown that, in this case, it is possible (although unlikely) 
to get more than one overflow per stage or a single overflow of more than one bit. 

It is impossible to get an overflow of more than two bits. (This can be seen from 
equations (57) since there are only two additions). 
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11.4.2 An Accuracy Analysis of the Fixed Point Calculation 


II .4.2.1 Introduction 

We will assume, in this analysis, that the inputs (i.e., the real and imaginary 
parts of X(j) or A(m)) are represented by B bits plus a sign. We assume the binary 
point lies to the left of the leftmost bit. We showed earlier that the magnitudes of these 
members of the array would generally increase as we moved from stage to stage in the 
calculation. Hence, the method of operation is to test for overflow within the inner 
loop. If there is no overflow, the calculation proceeds as usual. If there is an over¬ 
flow, then the two inputs producing the overflow are shifted right until there is no 
overflow. The amount of the shift is recorded (it will be either one or two bits) 
and the entire array is shifted right this same amount. In this scheme, wc shift 
not only those elements we've already calculated but also those yet to be done. 

(This latter shifting could be done in subsequent calculations if this were more 
efficient). The total number of shifts is accumulated and the power of two, raised 
to the negative of this total number of shifts, constitutes an overall scale factor 
to be applied to the final array. 

Our numbers in most applications will come from an A-D converter and will 
contain a quantization error. We will not consider the propagation of this quantiza¬ 
tion error. Our purpose is to compare the accuracy of a short word fixed point 
computation with a longer word floating point computation. In practice, the quanti¬ 
zation error would be common to both alternatives. 

There are two operations which produce errors which are propagated through 
the calculation. 

a. When two B bit numbers are multiplied together a 2B bit product results. 

If this product is rounded to B bits, an error whose variance is: 

a 2 _ 2~ 2B /12 (60) 

is created. This error has a standard deviation: 

A = 2" B //l2 ~0.3(2 _B ) (61) 

b. When two B bit numbers are added together and there is an overflow, 
then the sum must be shifted right and a bit lost. If this bit is a zero, 
there is no error. If it is a one, there is an error of +2"® depending 
upon whether the number is positive or negative. The variance of this 
error (it is unbiased assuming there are an equal number of positive and 
negative numbers) is: 

a 2 = 2 -2B /2> (62) 


60 


It has a standard deviation 


A 2 = 2~ B 1 «0.7(2 B ) 


(63) 


11.4.2.2 Upper Bound Analysis 


In this section, we give an upper bound analysis of the ratio of the RMS 
error to the R M S of the answer. This upper bound is obtained by assuming 
that during each step of the calculation there is an overflow and a need to rescale. 
We let X|< be a typical real element at the kth stage (i.e., the real or imaginary 
part of a complex element) and let 


V(XjJ = variance of X^. 


2 2 

We will, in this analysis, replace A^ by GA 


r 


We will also let 



(64) 


Since the first stage gives an overflow, the original data must be rescaled 
or truncated by one bit. Hence, 


V(X 0 ) = 6 A 2 . (65) 

In going from the original data to the results of the first stage, W=1 and, hence there 
is no multiplication and we either add or subtract. Further, we assume that the next 
stage will result in an overflow and hence we will have to rescale. This gives 


= 2V(X Q ) + 4 • 6 A 2 
V(X ] ) - 2(6A 2 ) + 4 • 6 A 2 . 


( 66 ) 


In these expressions and this entire discussion, we are assuming all errors to be 
independent and, hence, that the variance of the sum is the sum of the variances. 
Going from the first stage to the second stage, we have W = (-1) ■ ^ and again there 
are only additions and subtractions. Thus, with the rescaling 


V(X.,) - 2V(X 1 ) + 4 2 • 6 A 2 

- 2 2 (6A 2 ) + 2(4- 6 A 2 ) + 4 2 • 6 A 2 . 


(67) 


In going from the second stage to the third stage, we have multiplications and wc 
have them in all subsequent stages. In generating the third stage, half the inner 
loops have multiplications. Consider the first equation of equations ( 57 ). All the 
other equations are identical in terms of error propagation. Remember that X^ (i) 
is complex; 
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( 68 ) 


Re |X 3 (i)} = Re {X 2 (i)} + Re{X 2 (j)} Re{W} - Im{X 2 (j)} Im{W} 


Equation (68) yields, with rounding to B bits after the addition and with rescaling. 


V'(X 3 ) = V(X 2 ) + [Re 2 {X 2 (j)} + Im 2 {X 2 (j)}] V(W) 

+ [Re 2 (W) + Im 2 (W)] V(X 0 ) 


(69) 


+ (4 3 A 2 ) + 4 3 • 6A 2 

= V(X 2 ) + |X 2 (j)| 2 A 2 + V(X 2 ) + (4 3 A 2 ) + 4 3 • 6A 2 


In Equation (69), the first term is the variance of the first term of (68). The second 
and third terms of (69) are the variance of the full 2B bit products given by the 
second and third terms of (68). The fourth term of (69) is the result of rounding 
after the addition. The fifth term is the rescaling term. Finally, we saw in Equa¬ 
tion (58) that the average modulus squared of the complex numbers is increasing 
by a factor of 2 every stage. Hence, if we let K equal the average modulus squared 
of the initial array, i.e., 


N-l 



i=0 

then we have: 


V' (X 3 ) = 2V(X 2 ) + 2 2 KA 2 + (4 3 A 2 ) + 4 3 • 6 • A 2 


(70) 


Equation (70) would be correct for V(X^) if all the inner loops involved multiplica¬ 
tions. However, at this stage only half*of them do and hence: 



In the next stage, three quarters of the inner loops require multiplications and 
these multiplications get progressively more numerous as the stages increase. 
Hence, from here on, we will assume all stages have multiplications in all the inner 
loops. Thus, applying the above techniques, we get: 
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V(X 4 ) = 2 4 (6A 2 ) + 2 3 (6 • 4A 2 ) + 2 2 (6 • 4 2 A 2 ) + 2(6 • 4 3 A 2 ) + 6 • 4 4 A 2 
+ 2 2 KA 2 + 2 3 KA 2 + 4 3 A 2 + 4 4 A 2 ; 


( 72 ) 


and generally if M is the last stage: 

V(X M )= 2 M (6A 2 ) + 2 M_1 (6-4A 2 ) + . . .+ 2(6-4 M_1 a 2 ) 

+ 2 M+2 KA 2 + (M-3)2 M_1 KA 2 + 2 M_4 (4 3 A 2 ) 
+ 2 M_4 (4 4 A 2 ) + . . .+ (4 M a 2 ) 


= (1.5)2 M + 2 A 2 (l+2+. . .+2 M *) + (M-2.5)2 M *KA 


0 M+3 2 0 M+4 _ 0 M-4 2 

+ 2 A + 2 (1+2+. . .+2 ) a • 


(73) 


Or 


V(X M ) « (1.5)2 2M+2 A 2 + (M-2.5)2 M *KA 2 + 2 M+3 A 2 + 2 2M+1 a 2 
^2 2M+3 A 2 + (M-2.5)2 M - 1 KA 2 + 2 M+3 A 2 - 


(74) 


K is the average of the square of the absolute values of the initial complex array. 
Hence, applying Parseval’s theorem, Equation (55), the average of the square of 
the absolute values of the final array will be 2^K. What is most meaningful in 
this case, however, is the mean square of the real numbers which is 2^K/2. 
Hence, we have: 


V(X ) _ 2 M+3 A 2 (M-2.5) A 2 /2 2 3 A 2 

-— ~ - + - + - 

\T * 

2 m K/2 K/2 1/2 K/2 

and finally, 


RMS (error) 

RMS (result 


M+3 



K/2 


M+3 

2 2 2~ B (0.3) 

RMS (initial array) 


(75) 


(76) 
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Equation (76) gives an approximate upper bound for the ratio of the RMS 
of the error to R M.S. of the answer. Notice that this bound increases as the \ r N 
or 1/2 bit per stage. 


II.4.2.3 Lower Bound Analysis 

We will now obtain a lower bound for the ratio of R M S of the error to the 
R M S of the answer. We obtain this lower bound by assuming that there are no 
overflows in the calculation and, hence, no shifts of the array. In this case, 


V(X 0 ) = V(X 1 ) = V(X 2 ) = 0. (77) 

In the third stage, half of the inner loops involve a multiplication and hence: 

V(X 3 ) = (1/2) (2 2 KA 2 ) + 1/2 (A 2 ). (78) 

This can be seen by considering the first term of equations (69). The first term 
of (78) comes from the second term of the first of equations (69). The second term 
of (78) is caused by the rounding to B bits. Now, as before: 


V(X 4 ) = 2V(X 3 ) + 2 3 KA 2 + A 2 

= 2 2 KA 2 + 2 3 KA 2 + A 2 + A 2 . 


(79) 


Finally, 


V(X M ) . 2 M - 2 KA 2 ♦ (M-3)2 M - 1 KA 2 . 2 M -V + 2 M ‘ 5 A 2 

_M-b .2 .2 

+ 2 A + . . . + A 

= (M-2.5)2 M_1 KA 2 + 2 M_3 A 2 + (1+. . ,+2 M_5 ) A 2 


( 80 ) 


~ /*/r o C w, M_1 TrA 2 0 M-3.2 0 M-4 .2 

~ (M-2.5)2 KA +2 A + 2 A . 

As in subsection 2.2, the mean square of the array of resulting real numbers is 
2 M . k/ 2. Hence, we have: 


V * X M> 

2 M K/2 


(M-2.5) a 


A^/S 

K/2 


2 /c 
A / 6 

K/2 


(81) 


(M-2.5) A . 
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and 


RMS (error) 

RMS (result) 


(82) 


(M-2.5) 


1/2 


The lower bound increases as M = log 2 N. This is the rate of increase which 
has been observed for the floating point calculation. ’ 1 


II.4.2.4 Some Experimental Results 

An IBM 7094 program was written to perform a fixed point calculation using 
the fast Fourier transform algorithm, as described above. The program was 
capable of simulating a fixed point machine of any word size up to 35 bits plus a 
sign. Experiments were run with fixed point numbers of 17 bits plus a sign. This 
corresponds to B = 17 in the analysis of subsections II.4.2.2 and H.4.2.3. 

The experiments were performed in the following way. Floating point input 
was fixed to 17 bits plus a sign. This fixed input was then transformed with the 
fixed point program. The fixed point output was then floated. Next, the fixed 
point 17-bit input was floated and a floating point transform taken. Since this 
floating point transform uses a floating point word with a 27-bit mantissa, it was 
considered the correct answer. Finally, the R M S of the difference between the 
fixed point and floating point answers was taken. We also obtained the maximum 
absolute error and average error. 

Figure 27 contains the result ol transforming random numbers which lie between 
zero and one (placed in both the real and imaginary parts). In this and subsequent 
tests, three runs were made for every power of 2 from 8 to 2048. Since these ran¬ 
dom numbers have a DC component of one-half, the fixed point program must 
rescale at least every stage but one. Hence, one would expect the error to lie close 
to the theoretical upper bound as given by Equation (76). This theoretical upper 
bound is also plotted in Figure 27 and the results are seen to lie slightly above it. 

The R M S of the original array, K/2, is approximately 0.6. 

Figure 28 contains the results of transforming three sine waves plus random 
numbers between zero and one-half in the real part and all zeros in the imaginary 
part. Specifically: 

Re{X(j)} - 1/2[Y(j) + (1/2) sin (27r8j/N) 

+ (1/4) sin (27r4j/N) + (1/4) sin (27T8j/N)J 

Im{X(j)} - 0 
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Figure 27. Experimental Error Results: Random Numbers 
Between 0 and 1, B=17 
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Figure 28. Experimental Error Results: Random Numbers Plus 3 Sine Waves 
0 < Random Numbers < 1/2, B = 17 

























where the Y(j) are random numbers between zero and one. Again, there is a 
DC component of magnitude one-fourth and the array must be rescaled at least 
at every stage but two. Thus, one would expect these results to be lower relative 
to the theoretical upper bound than the case depicted in Figure 27. From Figure 28 
one can see that this is in fact the case. The R M S of the original array, K/2, 
is, in this case, approximately 0.35. This is the reason the upper bound curve is 
higher than that of Figure 27. 

Figure 29 contains the results of transforming random numbers from minus 
one to one (in both real and imaginary parts). In this case, the DC component is 
zero and there is no other strong component. The number of shifts should be 
approximately (log2N)/2 or one-half shift per stage. Hence, one would expect 
the error curve to lie well below the theoretical upper bound as is the case. In 
this case, K/2 = 0.6. 

Figure 30 contains the results of an experiment identical to that corresponding 
to Figure 28, except the random numbers are between plus and minus one half. The 
results are as expected. In this case, K/2 = 0.35. 

Finally, Figure 31 contains the results of transforming a sine wave in the real 
part and zero in the imaginary part. The sine wave was sin (27rj/8). Although in 
this case the array must be rescaled in at least every stage but two, the error is 
well below the upper bound. Here, K/2 = 0.5. 

These calculations were made with rounding rather than truncation. The error 
with truncation was somewhat larger. Rounding can be included in the microcode 
described in Section II.4.3 with a negligible decrease in speed. The bias, as 
reflected by the average error, was in general negligible compared with the 
RMS error. The maximum error was roughly what would be expected. 


II.4.2.5 Conclusions and Additional Comments 

Figure 32 below gives the upper bound on the ratio of R M S of the error to the 
R M S of the answer assuming B = 15 and assumming that the RMS of the initial 
array is three-tenths (i.e., K/2 = 0.3). This corresponds to the microcoding 

situation for the Model 40. 

There are two main uses to which the fast Fourier transform would be put in 
seismic signal processing. The first would be in estimating the spectral matrix 
of a subarray of seismometer outputs as a first step in the calculation of the coef¬ 
ficients for optimum processing (see 13, 14, 2 or 8). The second would be in 
estimating a sequence of short term spectra of the LASA beam output as one of 
the first steps in discrimination processing. In both of these cases, an acceptable 
and computationally efficient procedure is to take transforms of short, overlapping 
segments of the time series; calculate periodograms or cross-periodograms and 
time average over these to increase the statistical accuracy. (See Section 8 of 
ref. 11, and ref. 13.) 
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Figure 29. Experimental Error Results: Random Numbers 
Between -1 and 1, B=17 
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Figure 30. Experimental Error Results: Random Numbers Plus 3 Sine Waves 
- 1/2 < Random Numbers < 1/2, B = 17 
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Figure 31. Experimental Error Results: Single Sine Wave, 
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Figure 32, Error Upper Bound for 15-Bit Calculation 


















As concerns the fixed point accuracy, these remarks are important in two 
respects. First, since short segments are transformed, the fixed point error 
would be kept small. In the spectral matrix estimation, the size of the segments 
would be of the order of the number of coefficients in the optimum processing 
procedure. In the estimation of LASA beam spectra, it would be of the order of 
the reciprocal of the resolution desired (in terms of the Nyquist frequency as 
unity). In both of these cases, 256 would be a high figure for the length of the 
segments to be transformed. For segments of length 256, the fixed point error 
would be small compared to the statistical error in the spectral estimates. Its 
only effect might be a slight one on the available dynamic range. Second, the 
time averaging would reduce this fixed point calculation error as well as the 
statistical error. Hence, for the purposes of spectral estimation important in 
LASA signal processing, a fixed point calculation appears to be more than ade¬ 
quate. 

As a final comment, since the magnitude of the error has a sensitive and direct 
relationship to the number of times the array is scaled, it would be best to remove 
any substantial DC component which happened to be present. 
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II.4.3 Speed Advantage of the Microcoded FFT 


II.4.3.1 The Inner Loops 

M 

The FFT algorithm operates on an array of 2 complex numbers, taking M 
stages on the array. Each stage comprises 2^-1 passes through the inner loop 
computations, where each pass processes two complex points of the array. In 
terms of these points, called OPerand X and OPerand Y, and the complex root 
of unity U, the inner loop is either a Type I or Type II: 


Type 1: 


Complex 

OPX 

OPY 

= OPX + OPY 
= OPX - OPY 

Real 

ROPX 

ROPY 

= ROPX + ROPY 
= ROPX - ROPY 

Imaginary 

IOPX 

IOPY 

= IOPX + IOPY 
= IOPX - IOPY 

Type II: 



Complex 

OPX 

OPY 

= OPX + OPY• U 
= OPX - OPY•U 

Real 

ROPX 

ROPY 

= ROPX + (ROPY - RU - IOPY- IU) 

= ROPX - (ROPY - RU - IOPY- IU) 

Imaginary 

IOPX 

IOPY 

= IOPX + (ROPY - IU + IOPY- RU) 

= IOPX - (ROPY- IU + IOPY- RU) 

The Type I inner loop is used for all passes in the first and second stages 
through the array. Thereafter, its application is halved for each subsequent 
stage (see Table 2). 

The number of Type II inner loop executions in the total FFT process does 
not exceed the Type I inner loops unless M > 5. Actually, the ratio Type 11/ 
Type I is asymptotic to M/2, but this is not a good measure of the relative fre¬ 
quencies of the types unless M is large. Nevertheless, the Type II inner loop 
is the dominant factor in the FFT process, because it involves four real multi¬ 
plications and six real additions. Type I involves no multiplications, and only 
four additions. For this reason, our estimate has centered around the Type II 
inner loop. 
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Table 2 
INNER LOOP 


Stage 

Number of Type I 

Inner Loops 

Number of Type II 

Inner Loops 

1 

2 M-1 

0 

2 

2 M- 1 

0 

3 

2 M- 2 

2 M-2 

M-l 

22 

0 2 3 . „M-2 

M 

2 

2 +2 2 + . . . + 2 M_2 


II.4.3.2 Type II Inner Loop Process and Flow Chart 

M 

To maximize and preserve significance in the fixed point operands, the 2" 
complex numbers are normalized in the sense that: 

a. The high order bit (bit position 0) is a sign bit 

b. The binary point is assumed to be between bit positions 0 and 1 

c. At least one real or imaginary component of one of the complex 
points has the value 1 in bit position 1. 

With this representation of the complex points, the execution of inner loop 
passes may give rise to an overflow of one or two bit positions. Actually, a two 
position overflow has very low probability; and a three position or more overflow 
is impossible. 

M-l 

Consider the 2 inner loop passes in any stage. If an overflow of one bit 
position occurs first in pass K, this overflow divides the stage into two parts: 

a. Previously executed passes; that is, passes 1 to K-l 

M-l 

b. Subsequently executed passes; that is, passes K+l to 2 

Similarly, previous and subsequent passes may be defined with respect to a 
two position overflow. This will divide the stage into three parts, if the two- 
overflow initially occurs in a pass subsequent to the one-overflow. 

When an overflow of one or two positions occurs, the entire array must be 
shifted right one or two bit positions to preserve significance. This is most 
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expeditiously done by shifting the results of each subsequently executed pass in 
the stage before these results are returned to main storage; and, after the stage 
is completed, by returning to the beginning and shifting the results of all pre¬ 
viously executed passes. The latter process, of shifting previously executed 
results, will be called post shifting. An array is post shifted until an array ele¬ 
ment with a stop-address is detected. This is the address of an inner loop 
operand in the pass which first gave rise to the overflow. Saving these addresses 
is a bookkeeping operation which must be performed when an initial overflow of 
one or two places occurs. 

From the above, we see that the nature of the inner loop process depends on 
whether an overflow condition has occurred in a previously executed pass, and 
on whether a new or additional shift condition has occurred in the currently exe¬ 
cuted inner loop pass. All of the possibilities are given in Table 3. 

For each case in Table 3, we give a probability. These assigned probabili¬ 
ties are very rough estimates. They are not to be taken literally, except as an 
indication that, for purposes of estimating microprogram execution time, it 
suffices to consider only Cases A, D, and E. The effect of the low probability 
Cases B, C, F, and G is negligible. These conclusions are consistent with error 
analysis results. 


Table 3 

INNER LOOP OVERFLOW CONDITIONS 


Case 

Ovfl Positions 
Previous Pass 

Ovfl Positions 
Current Pass 

Probability 

Total Shift 
Before Storing 

A 

0 

0 

0.500 

0 

B 

0 

1 

0.010 

1 

C 

0 

2 

0.001 

2 

D 

1 

0 

0.230 

1 

E 

1 

1 

0.230 

1 

F 

1 

2 

0.001 

2 

G 

2 

0 

0.028 

2 


If two 16-bit normalized binary numbers (binary points assumed between bit 
positions 0 and 1) are multiplied together, as with SUMP hardware, the high half¬ 
word of the result will have its binary point between bit positions 1 and 2. Hence, 
the high halfword of the product cannot be added to a normalized number without 
first either shifting the product left one position, or shifting the normalized number 
right one position, to line up the binary points. The alternative chosen depends on 
the part of the inner loop being processed (computation of RX and RY, or computation 
of IX and IY), and on whether or not there is a pending shift from a previous pass. 
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A flowchart of the inner loop process, showing the testing and handling of 
the various overflow conditions, is shown in figures 33, 34, and 35. The states 
Y2 and Y3 are used to signal overflow conditions from previously executed 
passes as illustrated in Table 4. 


Table 4 

OVERFLOW CONDITIONS 


Y2 

Y3 


0 

0 

No Pending Overflow 

0 

1 

Pending Overflow of 1 Position 

1 

1 

Pending Overflow of 2 Positions 


In the currently executed pass, the occurrence of overflows in the computa¬ 
tion of RX and IX is recorded in the B1 REG, and overflows in the computation 
of IlY and IY are recorded in the BO REG. 

New overflows, conditions B, C, and F, can be recognized either in the cal¬ 
culation of RX and RY, or IX and IY. Thus, we have two exits each for these con¬ 
ditions in the flow charts. 

A three-character notation near the lower right of the flow chart boxes cor¬ 
relates the flow chart with the microprogram given on figures 36, 37, 38, and 39. 
The first character gives the sheet number, and the second two give the block 
coordinates of the microinstruction on the sheet. 


II.4.3.3 Execution Time Estimates 

The microcoded execution time estimates for the Type II inner loop cases 
shown in figures 34 and 35 are given in Table 5 and assume both fixed point 
operands and the SUMP hardware for multiplication. 

Figures 40 and 41 give the OS/360 program for the inner loops. The execu¬ 
tion for these is as follows: 

Fixed Point, assuming no pending shift: 485 psccs 

Fixed Point, assuming one pending shift: 480 psecs 

Floating Point: 560 psccs 
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Figure 33. Flow Chart for inner Loop of Fast Fouiier T-an form-f- 
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Concurrent Processing-*■ 


SUMP(RY • IU + IY • RU) 


STORE RX and RY in MS 

2J5 


2L1 


1 



(RY • IU + IY 

• RU) LEFT SHIFT 



3C5 



CASE A EXIT 3N3 


Figure 34. Flow Chart for Inner Loop of Fast Fourier 
Transform—II 
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CASE B 


CASES D, E 


Figure 35. Flow Chart for Inner Loop of Fast Fourier Transform—III 
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Figure 36a inner Loop Microprogram ! 
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F-gure 37„ Inner Loop M'c r oprog r am il 
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F:gu r e 38, nner Loop M c opiog'am HI 
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F>gure 39, Inner Loop Microprogram IV 




























Table 5 


TYPE II INNER LOOP EXECUTION ESTIMATES 


Case 

Probability 

Microinstructions 

Microseconds 

A 

0.5 


86 

54 

B 

0.01 

78 

+ (5)* 

52 

C 

0.001 

81 

+ (20)* 

63 

D 

0.23 


78 

49 

E 

0.23 


78 

49 

F 

0.001 

73 

+ (15)* 

55 

G 

0.028 

(85)* 

53 




Average 

52 


*In the cases not microcoded, the numbers in parentheses are estimates 
of the instructions to accomplish the following: 

B: Set Y3, save array stop-addresses for post shift processing 

C: Save stop addresses, shift operands 2 places, store 

F: Save stop addresses, shift operands 1 place, store 

G: Inner-loop process, shifting all results two places. 


85 








REGISTER 1 CONTAINS ADDRESS OF OPERAND X, AOPX 
REGISTER 2 CONTAINS ADDRESS OF OPERAND Y, AOPY 
REGISTER 3 CONTAINS REAL PART OF W,RW, LEFT ADJUSTED 
IN LOW HALF WORD 

REGISTER 4 CONTAINS IW, LEFT ADJUSTED IN LOW HALF WORD 
PROGRAM 1 FIXED POINT DATA 


LR 

5,3 

REG5=RW 

LR 

6,4 

REG6=IW 

LR 

7,3 


LR 

8,4 


MH 

5,0(2) 

REG5=RW :: RY 

MH 

6,2(2) 

REG6=IW”IY 

MH 

8,0(2) 

REG8=IW”RY 

TM 

PRIOR,X 1 

01' IF BIT 7 OF PRIOR IS 1 CONDITION CODE IS 1 

BC 

1, LET 

BRANCH TO LET IF AT LEAST 1 PENDING SHIFT 

NO PRIOR PENDING SHIFT 

SR 

5,6 

RW :: RY-IW :: I Y 

SLA 

5,1 


BC 

l,OVFL 

OVERFLOW CASE NOT PROGRAMMED 

AR 

7,8 

RW :: I Y+IW :: RY 

SLA 

7,1 


BC 

1,OVFL 

OVERFLOW CASE NOT PROGRAMMED 

SRA 

5,16 


SRA 

7,16 


LH 

9,0(1) 

REG9-OPERAND X REAL PART 

LR 

10,9 

REGIO^OPERAND X 

SR 

9,5 

RY=RX-(RW :: RY-I W :: I Y) 

BC 

l,OVFL 

OVERFLOW CASE NOT PROGRAMMED 

AR 

10,5 

RX=RX+(RW :: RY-1W«IY ) 

BC 

l,OVFL 

OVERFLOW CASE NOT PROGRAMMED 

STH 

9,0(2) 

STORE RY 

STH 

10,0(1) 

STORE RX 

LH 

9,2(1) 

REG9=IMAGINARY PART OF OPAND X 

LR 

10,9 


SR 

9,7 

I Y=IX-(RW"IY+IW :: RY) 

BC 

l,OVFL 

OVERFLOW CASE NOT PROGRAMMED 

AR 

10,7 

I X= IX+ (RW :: IY+1 W :: RY ) 

BC 

l,OVFL 

OVERFLOW CASE NOT PROGRAMMED 

STH 

9,2(2) 

STORE IY 

STH 

10,2(1) 

STORE IX 


END OF NO PRIOR PENDING SHIFT CASE 


Figure 40. OS/360 Programs and Execution Estimates for 
FFT Inner Loop (I, II, and III) 
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AT LEAST ONE PRIOR SHIFT PENDING 


LET 


SR 

5,6 

RW :: RY-IW :: I Y 


AR 

7,8 

RW :: I Y+IW :: RY 


TM 

PRIOR,X' 

02' IF BIT 6 OF PRIOR 

IS 1 CONDITION CODE 

BC 

1,ANTHER 

2 PRIOR SHIFTS,CASE NOT PROGRAMMED 



ONE PRIOR SHIFT 


SRA 

5,16 



SRA 

7,16 



LH 

9,0(1) 

REG9=RX 


SRA 

9,1 



LR 

10,9 

REG10=RX 


SR 

9,5 

RY=RX-(RW :: RY-IW :: IY) 


BC 

l,OVFL 

OVERFLOW CASE NOT 

PROGRAMMED 

AR 

10,5 

RX=RX+(RW :: RY-1 W :: IY ) 


BC 

l,OVFL 

OVERFLOW CASE NOT 

PROGRAMMED 

STH 

9,0(2) 

STORE RY 


STH 

10,0(1) 

STORE RX 


LH 

9,2(1) 

REG9=IX 


SRA 

9,1 



LR 

10,9 

REG10=IX 


SR 

9,7 

IY=IX-(RW :: IY+IW :: RY) 


BC 

l,OVFL 

OVERFLOW CASE NOT 

PROGRAMMED 

AR 

10,7 

IX=IX+(RW"IY+IW :: RY ) 


BC 

1,OVFL 

OVERFLOW CASE NOT 

PROGRAMMED 

STH 

9,2(2) 

STORE IY 


STH 

10,2(1) 

.STORE IX 



END OF ONE PRIOR SHIFT CASE 


Figure 40. (continued) 
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TIME FOR INNER LOOP FIXED POINT 
ASSUMING NO PENDING SHIFT 
INSTRUCTION NUMBER TIME TOTAL TIMES 


LR 


6 

7.3 

45.0 

MH 


4 

43.0 

180.0 

TM 


2 

9.38 

18.76 

BC 


6 

9.38 

56.28 

SR 


3 

7.5 

22.5 

AR 


3 

7.5 

22.5 

SRA 

16 

2 

18.75 

37.50 

SLA 

1 

4 

10. 

40. 

LH 


2 

18.13 

36.26 

STH 


2 

10.63 

21.26 

TOTALS 


34 


480.06 MICROSECONDS 


TIME FOR INNER LOOP FIXED POINT 
ASSUMING ONE PENDING SHIFT 


INSTRUCTION 

NUMBER 

TIME 

TOTAL TIMES 

LR 

6 

7.5 

45.0 

MH 

4 

45.0 

180.0 

TM 

2 

9.38 

18.76 

BC 

6 

9.38 

56.28 

SR 

3 

7.5 

22.5 

AR 

3 

7.5 

22.5 

SRA 16 

2 

18.75 

37.50 

STH 

4 

10. 

40. 

SRA 1 

2 

18.13 

36.26 

LH 

2 

10.63 

21.26 

TOTALS 

34 


480.06 MICROSECONDS 


Figure 40. (concluded) 
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REG1 CONTAINS AOPX,IMAGINARY PART AT AOPX+4 
REG2 CONTAINS AOPY,IMAGINARY PART AT AOPY+4 
FLREGO CONTAINS RW 
FLREG2 CONTAINS IW 


LER 

4,0 


FLREG 4 = RW 

LER 

6,2 


FLREG 6 = IW 

ME 

4,0(2) 


RW”RY 

ME 

6,4(2) 


IW”IY 

SER 

4,6 


RW"RY-IW”IY 

LE 

6,0(1) 


RX 

SER 

6,4 


RX-(RW :: RY-IW”IY) 

STE 

6, TEMP 


TY TEMP STORE 

LE 

6,0(1) 


RX 

AER 

6,4 


RX+(RW”RY-IW”IY) 

STE 

6,0(1) 


RX STORED 

LER 

4,0 


RW 

LER 

6,2 


IW 

ME 

4,4(2) 


RW :: I Y 

ME 

6,0(2) 


IW-RY 

AER 

4,6 


RW”IY + IW”RY 

LE 

6,4(1) 


IX 

SER 

6,4 


IX-(RW”IY+IW”RY) 

STE 

4,4(2) 


STORE IY 

LE 

6,4(1) 


IX 

AER 

6,4 


IX + (RW”IY+IW”RY) 

STE 

6,4(1) 


STORE IX 

MVC 

4(4,1), 

TEMP 

STORE RY 


END OF 

FLOATING POINT PROGRAM 

TIME FOR 

INNER LOOP 

FLOATING POINT 

INSTRUCTION NUMBER 

TIME 

TOTAL TIMES 

LFR 

4 

7.5 

30.0 

ME 

4 

80.63 322.52 

SER 

3 

14.3 

42.9 

LE 

4 

11.88 45.76 

STE 

4 

12.5 

50.0 

AER 

3 

14.3 

42.9 

MVC 

1 

26.25 26.25 

TOTALS 

23 


560.33 MICROSECONDS 


Figure 41. Floating Point Data 
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These figures, divided by the 52 /isecs microcode estimate, give: 


Advantage over OS/360 Fixed Point: 9.2 

Advantage over OS/360 Floating Point: 10.8 


Remarks on Above Results 

a. Note that the inner loop with 0-pending and 0-found overflow is longer 
in S/360 execution than the case with 1-pending and 0-found. This 

is so because the normalized format of the data gives a natural right 
shift of one place in multiplication. This must be unshifted in the 
no-pending shift case; in the 1-pending shift case the result is cor¬ 
rectly positioned after the multiply, without further processing. 

Taking advantage of this gives rise to the discrepancy between Case A, 
and cases D and E, making the pending shift cases D and E faster. 

b. The floating point S/360 process does not appear unduly longer than 
the fixed point (560 to 480). Note that the floating point process 
involves no testing for overflow, and no special provisions for shifting. 

c. In the microcoded cases a pending shift is processed by shifting the 
result of the innerloop computation. Alternatively, a pending shift may 
be executed on the input data before the inner loop processing. 

If input data is shifted, there is an improvement in speed over the output 
result shifting because some ALU idle time during SUMP operation 
can be utilized for the process. But the preshift loses accuracy: 

Note that if ct denotes shift, 


ctx + ay = <r(x+y) 

For o (01) + a (01) = 0 + 0 = 0 
but ct (01 + 01) = o-(10) = 1 

The amount of time saving in executing the pending shift prior to inner 
loop execution will be about five microinstructions, or about five per¬ 
cent of the inner loop. 

d. The estimates given are for the power-of-two FFT process. 
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11.4.3.4 The Overall FFT Process 


Besides the inner loop, the Fast Fourier Transform comprises: 

a. Index Processing, for accessing operands for the inner loop 

b. Initialization, instruction fetch and setting initial counter values 
and working local storage. 

c. Post shifting, suppose an operand array has 256 data points. The 
inner loop is executed 128 times in each stage. If an overflow occurs 
first in the 64th execution, all inner loop outputs from the 64th to the 
128 will be properly positioned in storage. But the stage will not be 
completed until the outputs from the first 63 inner loop executions 
are also right shifted, in the post shifting process. 

d. Data scramble routine, for bit reversal of operand addresses. A 
FORTRAN program (Figure 42), flowchart (Figure 43) and microcode 
estimate (Figure 44) of the scramble routine are given. The micro¬ 
programming is illustrated by figures 45 and 46. It is interesting to 
note that the advantage for the scramble routine is appreciably higher 
than the usual seven or eight for a scientific application. In general, 
the advantage of microcode derives from the elimination of I-Fetches 
in a sequence of system instructions, and from the determinate oper¬ 
ands. In the microcoded scramble routine, there is another factor: 

The "branch on condition" is virtually free in microcode, but costs a 
system instruction. The scramble routine's high advantage results 
from the highly branching nature of the routine. 

It is well known that the number of inner loop executions in the total FFT 
processing of 2^ = N complex points is M • 2^ - l (M is the number of stages 
through the array, and 2 points of an array are processed with each inner loop 
execution). Experience indicates that the overall FFT execution time may be 
estimated as twice the inner loops execution time: 

2 • M • 2 M_1 - t • N • log 0 N, 

where the initial 2 provides for items a through d above, and where t is execu¬ 
tion time per inner loop. Our estimate of t is 52 microseconds, so that: 


52N log^N, usees 

is the microcoded FFT process time for an array of N points. A few values are 
given in Table 6. 
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J=1 

Nl=N/2 

IMAX=N-1 

DO 3 1=1,IMAX 

IF(I.GE.d) GO TO 3 

XT=X(J) 

X(J)=X(I) 

XI=XT 

YT=Y(J) 

Y(J)=Y(I) 

Y CI)=YT 

5 K=NT 

4 IF(K.GE.J) GO TO 3 

J=J-K 
K=K/2 
GO TO 4 

3 J=J+K 

RETURN 


{ 


OS/360 PROGRAM 


INIT 

LH 

2, ONE 

REG2=J,INITIALLY = 1 



LH 

3,N 





SRA 

3,1 

REG3=NT, INITIALLY = 

N/2 


LH 

8, ONE 

LOAD INCREMENT 



LH 

1, ONE 

REG1=I, INITIALLY = 

1 


LH 

9,IMAX 

LOAD COMPARAND 



LR 





LOOP1 

SR 

4,2 ' 

l\LOOP 

2 I-J 



BC 

4,EACH 

/ 



LOOPX 

LR 

5,3 

\ LOOP 

1 K=NT 


LOOP3 

r LR 

6,5 | 

f 




SR 

6,2 


K-J 


) 

BC 

10,UPIjJ 

1 



< 

SR 

2,5 


J=J-K 



SRA 

5,1 





l BC 

AR 

13,LOOP3 



UPIJ 

2,5 


J=J+K LOOP 

1 


BXLE 

1,8,LOOP1 

1 = 1+1 



EXIT 





EXCH 

LR 

11,1 

> 




SLA 

11,2 





AR 

11,BASE 





LM 

13,13,0(11) 

\ LOOP 2 



LR 

12,2 





S LA 

12,2 





AR 

12,BASE 





MVC 

0(4,11), 

,0(12) 

SET NEW XT 



STM 

13,13,0(12) 

SET NEW YT 



BC 

15,LOOPX j 





Figure 42. Jennrich Scramble Fortran Program 
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Loop 1 


N = 2 M 


A(0) = MS (BASE + 4) 


Initialization: 

j = i 


NT = N/2 


1 MAX = N 


BASE = BASE - 4 




1 = 

= 1 





Q-GzD-i O 



K = 

= NT 





gm k ~ . j y ~o 



J = J + K 


J = J - K 

1 = 1 + 1 


K = K/2 


I 



EXIT 


XT = BASE + 41 
YT = BASE + 4J 
Read Write MS (XT) 
BREG = DREG 
Read = MS(YT) 
CREG = DREG 
DREG = BREG 
Write MS(YT) 

Read MS(XT) 

DREG = CREG 
Write MS(XT) 

. I - 


Figure 43. Jennrich Scramble Flow Chart 
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LOOP 1 

I-J GREATER THAN OR EQUAL TO 0 
K-J GREATER THAN OR EQUAL TO 0 
INSTRUCTION NUMBER TIME TOTAL TIMES 


SR, AR 

3 

7.3 

22.5 

BC 

2 

9.38 

18.76 

LR 

2 

7.3 

15.0 

BXLE 

1 

16.26 

16.26 

TOTALS 

8 


72.52 MICROSECONDS 

LOOP 2 
INSTRUCTION 

NUMBER 

TIME 

TOTAL TIMES 

AR, SR 

3 

7.5 

22.5 

BC 

2 

9.38 

18.76 

LR 

2 

7.5 

15. 

SLA 

2 

25. 

50. 

LM 

1 

6.25+5 

11.25 

MVC 

1 

16.25+2 

. 5 :: 4 26.25 

STM 

1 

6.87+5 

11.87 

TOTALS 

12 


155.63 MICROSECONDS 

LOOP 3 
INSTRUCTION 

NUMBER 

TIME 

TOTAL TIMES 

LR 

1 

7.5 

7.5 

SR 

2 

7.5 

15. 

BC 

2 

9.38 

18.76 

SRA 

1 

18.13 

18.13 

TOTALS 

6 


59.39 MICROSECONDS 


SUMMARY 


LOOP 

S/360 CODE 

MICROCODE 

MICROCODE TIME 



MICROSECS 

INSTRUCTIONS 

MICROSECS 

ADVATAGE 

1 

72.52 

10.5 

6.56 

11 

2 

155.63 

22 

13.75 

11 

3 

59.39 

6 

3.75 

15 


Figure 44. Speed Advantage Estimates for the Loops 
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INTERFACE 

IN LOCAL STORAGE: J*!,NT*N/2 > I MAX 
BREG*I DBEO-J 


u. ^ 



K 

m x 


r 
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Figure 45, Jennrich Scramble Microprogram I 
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Figure 46. Jennrich Scramble M'croprogram II 
















Table 6 


FFT TOTAL TIME 


Array Size 

Total Time 

M 

N 

4 

16 

3 m secs 

8 

256 

100 m secs 

12 

4K 

2.3 sec 

16 

64K 

52 sec 


II.5 ARRAY OPTIMIZATION 


II. 5.1 Introduction 

An aspect of array design based upon maximizing gain as a result of judicious 
choice of element placement in accordance with a model of the noise field is pre¬ 
sented. The technique has utility in a relative sense since its practical application 
is dependent upon and limited to situations where noise data is available and a 
comparative analysis is desired. In Section 2, the formulation and procedure is 
outlined and some results are presented in Section 3. 


II. 5.2 Formulation 

Assuming signal coherence and equal noise power at each element in an N 
element array, the conventional processing gain*’ can be expressed as: 


G 


10 log 


N 


10 


1 + (N-l) p 


-dB. 


(83) 


where the average correlation coefficient is: 


P 


N(N-l) 


N 

I 

i,J=l 


P ij 


(84) 


and where p^ 


is the pair-by-pair element noise correlation coefficient. 
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By replacing with p(rj.), data derived models of the correlation coeffi¬ 
cient as a function of distance Between elements (r^) can be synthesized for 
application in Equation (83). By restricting the class of functions to prohibit 
maximums, the average value of the differential of p is: 


N 

d^ =(i/n)^s. • dr. , 

i=l 

where 


2 

N 

y 

~32<r >~ 


r.. 

ij 

N(N-l) 

L 

j=i 

J 


L 1 ij U 




(85) 


( 86 ) 


The arrow notation is used to denote vector quantities. 

A 2 

To simplify computation p has been modified to be a function of r.. as 

follows: ^ 


S. 

i 


N a, 2, 

4_ V 8 P (r ij ) - 

N(N-l) l 2 r ij 

j=l ij 

jA 


(87) 


This eliminates the need for computing | 7jj |. If "ry equals zero when computing 
tq, instrument j is moved to reduce the magnitude of the rj which is measured 
relative to the center of the array. This restricts two instruments from having 
the same position.^The maximum change in p is produced when drq is in the 
same direction as Sj; therefore, the sensitivity vector Si computed for each 
instrument determines the direction of movement. The instrument movement is 
determined by: 
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where K controls the incremental step size, and, thus, the magnitude of each 
instrument move is proportional to the ratio of its sensitivity vector to the sys¬ 
tem sensitivity. 


II.5.2.1 Constraints 

If an instrument is on the array constraint boundary, an instrument movement 
outside of the boundary is modified to be the projection of the Aq on the tangent 
to the constraint circle at the instrument position. When an instrument is moved 
outside of the constraint boundary, the instrument is returned to the constraint 
boundary along the instrument vector relative to the center of the array. A simi¬ 
lar technique is used for other types of constraints. For a spoke configuration 
where all instruments must lay on spokes through the center of the array, in addi¬ 
tion to being within a fixed diameter circle, the instrument position vector is 
rotated to correspond to the nearest spoke after each move. 


II.5.2.2 Program Operation 

After computing the directions of movement, all instruments are moved at 
the same time. Then, a new set of sensitivity vectors is computed. This process 
continues until either the sensitivity vector or its projection on the constraint 
(when constrained) falls below a variable threshold. A three instrument configu¬ 
ration is shown in Figure 47 to demonstrate the vector relationship associated 
with the program. For sufficiently small values of K, this technique should pro¬ 
duce convergence along the path of steepest descent toward an array configuration 
with a minimum value of p*. 


II.5.2.3 Operator Control 

To permit the array designer to dynamically control the threshold and incre¬ 
ment of instrument position change, the array configuration is displayed on an 
IBM 2250 CRT Display. Controls for modifying the program parameters are 
provided. 

The array configuration that satisfies the threshold setting is a function of 
the threshold values. A better instrument configuration can be achieved by 
reducing the magnitudes of both the threshold and instrument position change. 

The array configuration thus achieved satisfies a local minimum value of p*. The 
program does not have a test to determine that the minimum corresponds to a 
global minimum (the ’’optimum" configuration). Confidence in selecting an array 
configuration can only be achieved by repeated runs (preferably with different 
initial conditions) to test for better design configurations. This technique is 
intended only to provide insight into selecting an array geometry by providing a 
model to serve as a departure point. A final array design must be tempered with 
geographic constraints and data acquisition considerations. 
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II.5.3 Results 


II.5.3.1 Array Performance 

Array configurations generated using the optimization program are shown 
in Figure 48 for circular constraints with diameters of 7 and 20 kilometers using 
an exponential function for the estimated noise correlation model. Predicted array 
performance using exponential functions with different decay rates are shown in 
figures 49 and 50. The configurations generated by the array optimization pro¬ 
grams are compared with the current LASA subarray configurations and configu¬ 
rations where all instruments are constrained to be on a circle. For the LASA 
seven-kilometer subarray, the inner rings are removed to demonstrate the change 
in performance. Using both noise correlation models, the LASA E3 subarray con¬ 
figuration approached the maximum gain obtained with the array optimization pro¬ 
gram. The LASA subarray configurations are shown in figures 51 and 52. 

The 16 kilometer NORSAR subarray with a modified hexagonal geometry 
currently being installed is shown in Figure 53. Perturbations in the geometry 
resulted from practical instrument siting considerations. A modified circular 
geometry and the configuration generated by the optimization program with a 
16 kilometer constraint circle, are presented in figures 54 and 55, respectively. 

The modified circular geometry also resulted from practical considerations. In 
Figure 56, the optimal gain geometry is superimposed on the NORSAR subarray 
to demonstrate the differences in the configurations. Table 7 summarizes the array 
performance for the three configurations using an exponential noise correlation 
model. The NORSAR subarray and the circular configuration are evaluated over a 
range of exponential decay rates. In addition, the minimum communication line 
length required in each of the configurations is tabulated. 


II.5.3.2 Measured Noise Correlation Data 

To verify the validity of the noise correlation models, measurements of IASA 
subarray performance are presented in figures 57, 58, 59, and 60. The subarray 
performance for both unfiltered and filtered cases is shown. These performance 
characteristics we obtained by measuring the average correlation coefficient for 
the subarray over an 80 second interval. The two filters used were third order 
Butterworth with bandpasses of 0.9 to 1.4 Hz and 0.7 to 3.0 Hz, respectively. 

Figure 61 shows the exponential noise correlation models used for evaluating array 
performance. These models can be compared with the measured correlation coef¬ 
ficient between instruments as a function of the distance between them for three 
noise records. An 80 second noise record was used for each case. The separate 
noise records are identified by the event name or location which they preceded. 

The correlation coefficients were fitted with a sixth order polynomial to generate 
a continuous function. This data is presented in figures 62 through 70. 
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Figure 48. Array Configurations for 20 Instruments 
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Figure 49. Predicted Array Performance 


103 



















Noise Reduction dB 


12 


10 


- T - 

Predicted Array Performance 
Using A Correlation Function 
Estimated by the Exponential 

p = exp(-0.72 |r.. |) 








N 







D = 7 Km_ 










A. 

Array Optimi 
with 20 Km D 

nation Configurations 
diameter Circle Const 

ralnt . 



B. 

C. 

D. 

Montana LASA E-3 Subarray 

All Instruments on Circle with 20 Km 
Diameter 

All Instruments on Circle with Diameter 
of 7 Km Corresponded to Optimum 
Configuration 

Montanna LASA Subarray Configuration 
Adding Inner Rings 

1 _ 1 _ 1 _ 



E. 


4 6 

N(Number of Instruments) 


20 25 


Figure 50. Predicted Array Performance II 
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Figure 51. 25 Seismometer E3 LASA Subarray 
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Figure 52. LASA Subarray 
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Figure 53. "Hexagonal" Geometry 
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Figure 54. Modified Circular Geofnetry(D w 16 Km) 


108 






D = 16 Km 
K = 1.1 Km 


8 

6 k 


v 0 


-2 

-4 


-6 

-8 



-6 -4 -2 


0 

Km 


Figure 55. Optimum Gain Geometry 
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Figure 56. Superimposed Geometry (D ~ 16 Km) 
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Figure 57. Observed Subarray Performance, Montana LASA E3 
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Figure 58. Observed Subarray Performance, 
Montana LASA I, 7 Km Subarray 
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Figure 59. Observed Subarray Performance, 
Montana LASA JJ, 7 Km Subarray 
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Figure 60. Observed Subarray Performance, Montana LASA III, 7 Km Subarray 
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Figure 61. Noise Correlation Model 






















Noise Prio- to Kurile Event of 11/21/66 




o 

< ^r 

>*— 

o i 
t O. 
o 

_Q Q- 
D CO 


in 

«88 

£o 

o . 

II 


CL 0° 

- II H 


-O > S-? 
Oo^- 

XX 
~o o o 



NO 


in 


co 


o • ^ > 

5^-00 

-Q Q- _c -o 
Deo t i: 
l/) CQvQl/1 


CN 


F 


<U 

U 

C 

o 

O 


oo <) N; 

o o o 


CN O CN 

O O O 

1 


^uapijpo^ UO|4D|9JJO^) 


117 


Figure 62. Correlation Coefficient vs Distance Between Instruments I 















Noise Prior to Kurile Event of 11/21/66 
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Figure 63 Correlation Coefficient vs Distance Between Instruments II 
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Figure 64. Correlation Coefficient vs Distance Between instruments III 
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Figure 65. Correlation Coefficient vs Distance Between Instruments IV 
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Figure 66. Correlation Coefficient vs Distance Between Instruments V 
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Figure 67. Correlation Coefficient vs Distance Between Instruments VI 
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Figure 68. Correlation Coefficient vs Distance Between Instruments VII 
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Figure 69. Correlation Coefficient vs Distance Between Instruments VIII 
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Figure 70. Correlation Coefficient vs Distance Between Instruments IX 

















Appendix III 


PHASE DELAYS 


Operation of the LASA beam former requires the use of 21 time delays—one 
for each subarray. Previous studies have demonstrated that plane wave steering 
of subarray beams must be augmented by correction factors to take account of 
travel time anomalies. 

A key requirement for both Detection and Event Processing is the library 
organization of phase delays. The implications of using a library form of process 
structure are examined in this appendix. Procedures, based on correlating a 
LASA beam with each of the subarray beams used in forming it, have been 
developed for calculating improved steering delays, given an initial set. These 
procedures described in Appendix II—1 will be employed to maintain, correct, and 
enlarge the library of phase delay parameters. Regions and sectors in U-AZI-space 
are defined, and associated subarray correction factors are computed, reflecting 
the discrepancies between theoretical and actual travel times observed for a 
sample of several hundred events. A test performed to evaluate the effectiveness 
of the subarray correction factors in reducing mis-steering loss is also described. 

The methodology to cover an arbitrary location, U, in inverse velocity space 
by a judicious choice of preformed, presteered subarray beams is addressed in 
Section 3. Special attention is given to subarray beam assignment in order to 
minimize the expected loss. 


III.l PROCESS STRUCTURE 

The Detection Processor forms LASA beams by summing in phase appropriate 
beams chosen from a set of preformed, pre-steered subarray beams. Since the 
delays required for steering each subarray are a function of the location of the 
LASA beam being formed, they will be obtained from a library of phase delays, 
organized by LASA beam coordinates. In Event Processing the library delays 
will be employed or the optimal delays will be obtained from correlations that can 
be generated for events of special interest. 
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III. 1.1 Phase Delay Library Parameters 

The entries for the phase delay library parameters will be a combination of 
several source factors among which are 

a. Regular plane wave steering delays 

b. Sets of subarray steering corrections each of which is valid over a 
specified region of U-AZI space 

c. Inter-regional corrections calculated from the sets of subarray 
corrections 

d. Optimal delays obtained as a result of an iterative process based on 
correlation methods. 

A library structure demands that requirements be specified for establishing, 
maintaining, and updating library entries. 

It is planned to employ plane wave delays for beam steering, particularly for 
steering to aseismic locations. 

For seismic areas, it is planned to modify the plane wave delays by means 
of empirically derived subarray correction factors; each set is valid over a speci¬ 
fied region in U-AZI space. The development of 35 sets of regional correction 
factors will be described in this appendix. 

Formulas have been developed for interpolating between regions in order to 
obtain correction factors for events located outside regions. Plane wave delays 
will be replaced by optimal delays for certain regions and for weak events. These 
will be obtained by use of the correlation program described in Appendix II-1. 

It will be routine practice to check the subarray regional corrections and the 
interpolation formulas on a regular basis. This will detect any systematic devia¬ 
tions which may reflect trends in the correction factors. 

Analysis of strong events will be the chief instrument used to correct and 
enlarge library phase delays. Detection and identification of a strong event will 
be followed by the input of its delays to the correlation program, which will calcu¬ 
late optimal delays for the event. The optimal delays may be used to check 
existing region correction factors and associated interpolation formulas, or may 
serve as the initial set of correction factors for a new region. 

Criteria must be established for determining when the correlation program 
is to be used for updating library entries, and for determining when an event may 
be considered so strong that its delays may be used to correct library parameters. 
The interplay between library parameters and the correlation program is shown 
in Figure 71. For routine events, appropriate phase delays will be obtained from 
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Figure 71 . Interaction between Phase Delays and the 
Correlation Program 
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the library. For special events, in particular, weak ones, the correlation program 
will be employed to calculate optimal delays, using as input either modified plane 
wave delays or delays for the beam on which the event was detected. 


III.2 ANOMALY CHARACTERIZATION 


III. 2.1 Introduction 

This section describes the development of steering correction factors 
resulting from an examination of several hundred events, each of which yielded a 
set of 21 time delay corrections. It is shown that these sets of correction factors 
can be grouped to produce average sets of corrections which are valid over judi¬ 
ciously selected geographic regions. Because of the irregular space and time 
distribution of earthquakes, these regions do not cover the entire earth. Conse¬ 
quently, it is necessary to interpolate to obtain correction factors for beams that 
are being steered to areas located between the nominally selected regions. The 
development of such interpolation formulas is described in this appendix. 


III.2.2 Method of Correction 

Theoretical steering delays must be amended by adding some correction fac¬ 
tor. This is currently being accomplished by defining for each event region and 
receiver site a travel time anomaly. This anomaly is the difference between the 
actual travel time and a theoretical travel time, based on a worldwide travel time 
average for events at a given range, such as those found in the Jeffreys Bullen 
Seismological Tables. The required event range can be obtained for each event 
from an event location based on the world seismic network. Time of origin is 
also determined by this network. The average anomaly for each receiver for a 
given geographic region is then calculated from a number of events occurring in 
that region, and is used to provide a correction to the theoretical travel time 
tables to give the required beam steering delays. The LASA system steering 
error for a given event would depend on how well the event arrival time differences 
fit the region averages. An interesting and informative study of time anomalies 
using this method on the LASA array has been completed by E. F. Chiburis.^^ 

Another possible method of organizing the time delays is to establish a best 
fitting plane or second order wavefront for each event, and to use as anomalies 
the variabilities of the deviations from this best fitting wavefront. As in the 
currently used method, the delays are based on the anomaly averages for a given 
region. 

A program to implement the wavefront method has been written. This pro¬ 
gram, called the Least Squares Plane and Quadratic Wavefront Program, accepts 
as inputs the positions, arrival times, and computation weights of up to 25 
seismometers (or subarray beams), and calculates three different least squares 
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wavefronts with the actual arrival time deviations from the best fits. The 
three fits are 

a. A best-fit plane wavefront 

b. A best-fit plane wavefront followed by an adjustment of arrival 
times to consider the expected effect of range upon wavefront 
curvature, with the above process repeated for any desired 
number of iterations 

c. A best-fit quadratic wavefront plus the three derived second 
derivatives or curvatures, expressed in array coordinates, 
range coordinates, and the principal coordinates of the quadratic 
surface. 

The input data for the wavefront program consisted of optically read arrival 
time data, made available by Seismic Data Laboratories (SDL) of Teledyne Indus¬ 
tries, Inc., and the VELA Seismological Center (Air Force Technical Application 
Center), Alexandria, Virginia. A total of approximately four hundred event 
arrivals were made available in their original form by SDL.l7 


ni.2.3 Procedure for Obtaining Steering Delay Corrections 

Chiburis describes the grouping of several hundred events into 36 geographic 
regions. Table 8, which uses the regions so defined together with later data from 
SDL, lists the number of events in each region. As described in Reference 16, a 
region is formed by grouping teleseisms whose subarray travel time anomalies 
are consistent within the region. The 377 events and the 36 regions form the 
starting point of the development reported herein. 

The 377 events have been processed through the Least Squares Wavefront 
Program and the resulting computed values of azimuth (AZI) and inverse phase 
speed (U) have been plotted for each event as shown in Figure 72. The other out¬ 
put of the wavefront program—the deviation in seconds of the actual arrival from 
the best fitting plane wavefront—serves as input to the Seismic Steering Delay 
Anomalies Program. Given a group of events, each with its set of 21 correction 
factors or deviations (one for each subarray), this program calculates the average 
deviation at each subarray. The difference between each event f s deviation and the 
average deviation of the group at each subarray is computed, yielding a set of 21 
differences, i.e., deviations of deviations. The average and standard deviation of 
this set of differences and the mis-steering loss in dB for the event's wavefront 
are calculated. 

The division of the 377 events into 36 regions^ furnished the initial groups 
used in the Steering Delay Anomalies Program. The program results indicated 
that some of the events yielded excessive dB losses and it was decided to regroup 
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Table 8 


377 EVENTS GROUPED INTO 36 REGIONS 


Region 

No. 

Geographic No. of Events 

Area In Region 

Region 

No. 

Geographic No. of Events 

Area in Region 

1 

Central-So. 
Alaska 

6 

19 

N. Atlantic 

2 

2 

Kodiak Island 

11 

20 

Azores 

2 

3 

Alaska Pen. 

3 

21 

No. West Indies 

5 

4 

Unimak 

4 

22 

Venezuela-So. 

West Indies 

3 

5 

Fox-Andreanof 

19 

23 

Eastern Mexico 

9 

6 

Andreanof 

4 

24 

N. Centr. Ameria 

8 

7 

Andr-Rat 

22 

25 

S. Centr. America 


8 

Near Island 

14 

26 

N. Columbia 

12 

9 

Kamchatka 

12 

27 

Peru-Brazil 

8 

10 

Kurile 

29 

28 

Peru-Bolivia 

13 

11 

Sakhal in- 
Hokkaido 

8 

29 

N. Chile- 
Bolivia 

16 

12 

Honshu 

Marianas 

20 

30 

N. Chile- 
Argentina 

12 

13 

Marianas- 

Carolines 

13 

31 

C. Chile- 
Argentina 

13 

14 

Taiwan-Ryuku 

10 

32 

West Mexico 

4 

15 

Kazakh 

9 

33 

Easter Island 

3 

16 

Hindu Kush 

7 

34 

Easter Island 

7 

17 

East Caucasus 

3 

35 

Samoa-Tonga 

22 

18 

Greece-Turkey 

14 

36 

Fiji Islands 

23 


131 









Figore 72. 377 Wavefronts, Calculated U and AZi as Observed 
at LASA 
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the events into new regions. There are 35 regions in the new grouping. They 
are plotted in Figure 73 and their coordinates are listed in Table 9. The revised 
regions are quite similar to the original 36. For example, revised regions 109 
through 123 correspond exactly to original regions 12 through 26 respectively. 

After regrouping the 377 events into 35 new regions, averages were again 
formed and the anomalies program was used to check the assignment of the 
events to the regions and to assign weights to the events. For those events 
which continued to show a loss of more than 1 dB at 1 Hz, a weight of zero was 
assigned while the others were weighted one. New averages were formed for 
each region, using only those events weighted one, i.e., using those events whose 
loss was less thanldB at 1 Hz. The anomalies program was used to provide a 
final check. This time there were no further changes—the events weighted one 
continued to show a loss of less than 1 dB at 1 Hz while those weighted zero 
showed even greater losses (as was to be expected since they had been excluded 
from the average forming prooess the second time). Thus, the phase delay aver¬ 
ages for the revised regions are arrived at by a two step iterative process, based 
on stable, similar events each of which experiences a loss of less than 1 dB at 1 Hz 
due to steering to these regional averages. Table 10 lists the number of events 
with weight one for each region. As shown in the table, there is a total of 335 
weight one events out of a total of 377 events and only 42 events were weighted 
zero. 

The Anomalies Program yields the average loss (indB) per region, i.e., the 
standard deviations of the events in a region are averaged, and converted to a 
loss. 


For an event occurring in a region, this average loss figure represents the 
loss in dB which results from using the regional averages as steering correction 
factors rather than a set of corrections that would be optimum for the particular 
event. Table 11 lists the average loss per region for the 35 revised regions. As 
can be seen from the table, che average loss is less than 1 dB at 1 Hz for every 
region. These results indicate that the regional averages can be instrumental in 
providing correction factors for beams deployed within regions. 


in. 2.4 Definition of Sectors 

To obtain correction factors for areas having no history of teleseismic 
activity, some form of interpolation between regions must be developed. Inter¬ 
polation between regions is based on the assumption that correction factors are 
functions only of azimuth or inverse phase speed or both. This implication or 
assumption must be considered somewhat questionable, for examination of the 
regional corrections reveals no discernible pattern in the average deviations. 
They appear not to be simple functions of azimuth or of range, but they fluctuate 
in an erratic manner. Figures 74, 75, and 76 illustrate the deviations per region 
for three specific subarrays. 
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egic 

No. 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 


Table 10 


NUMBER OF CONTRIBUTING EVENTS 


No. of 
Events 

6 

10 

7 
37 
13 
16 
13 
16 
19 
11 

8 
9 


Region 

No. 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 


No. of 
Events 

5 

2 

12 

2 

2 

5 

3 
9 
8 
7 

10 

4 


Region 

No. 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 


No. of 
Events 

8 

6 

8 

13 

10 

7 

4 

5 
2 
4 

34 
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Table 11 


AVERAGE LOSS PER REGION 


Region 

No. of Events 

Weight One/Total 

Average Loss 
(in dB) 

101 

6/6 

0.44 

102 

10/13 

0.55 

103 

7/7 

0.48 

104 

37/41 

0.45 

105 

13/16 

0.34 

106 

16/17 

0.39 

107 

13/14 

0.33 

108 

16/18 

0.25 

109 

19/20 

0.36 

110 

11/13 

0.34 

111 

8/10 

0.31 

112 

9/9 

0.39 

113 

5/7 

0.49 

114 

2/3 

0.43 

115 

12/14 

0.45 

116 

2/2 

0.15 

117 

2/2 

0.15 

118 

5/5 

0.36 

119 

3/3 

0.77 

120 

9/9 

0.53 
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Table 11 


Region 


121 

122 

123 

124 

125 
12G 

127 

128 

129 

130 

131 

132 

133 

134 

135 

Total 


(concluded) 

No. of Events 
Weight One/Total 

8/8 

7/7 

10/12 

4/4 

8/9 

6/7 

8/8 

13/14 

10/13 

7/7 

4/4 

5/7 

2/3 

4/5 

34/40 

335/337 


Average Loss 
(in (IB) 


0.45 

0.27 


0.38 


0.36 


0.41 


0.36 

0.32 


0.39 


0.44 

0.41 

0.21 


0.54 


0.16 

0.41 


0.36 



Figure 74 Regional Deviations (seconds) for Subarray Cl 
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Figure 75. Regional Deviations (seconds) for Subarray C2 
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An interpolation procedure has been developed in an attempt to estimate 
correction factors for areas near regions. The interpolation method is based 
on the construction of areas called sectors which are essentially defined in terms 
of regions. A sector is an area containing two or more regions. The guideline 
used in establishing sectors has been to define their boundaries so that a sector 
is essentially constant either for azimuth or inverse phase speed. Figure 77 
depicts the nine sectors that have been defined. 

Table 12 gives the U-AZI coordinates tor each sector and lists the regions 
contained in each sector. As shown in Figure 78 sectors do not extend significantly 
beyond the regions they contain. 


III.2.5 Methods of Interpolation and Extrapolation 

The interpolatory method used to obtain a deviation for each subarray and 
sector is simple linear interpolation in which either azimuth or speed is the 
independent variable. For example, if an event occurs in sector A, the deviation 
to be used for steering to that event is obtained by interpolating in azimuth between 
regions 114 and 115. The linear interpolation formula is stated in Table 13. 

Table 14 gives the (X, Y) values relative to subarray B1 for each region in each 
sector. Here, Y is the value of the deviation while X is either an azimuth value or 
a U value, depending on whether the independent variable is azimuth or inverse 
phase speed. The tables for the other subarrays are similar and therefore not 
presented. 

Two or more regions are sometimes combined and treated as a single region. 
In particular, regions 122 and 124, which are contained in Sector C, are assigned 
the common azimuth value of 147° and the common deviation of 0.017 seconds; this 
pair of values is used for all 21 subarrays. 

The interpolation procedure will be demonstrated by computing the correction 
factor for subarray B1 steered to event V whose location is estimated at U = 0.080 
and AZI = 130°. Figure 78 shows that V is in sector C but not within any region 
in C. Table 14 shows that deviation is a function of azimuth in sector C; with 
respect to that variable, V is located between regions 119 and 123. Employing 
the notation of Table 13, let (X^, Y^) represent the coordinates of region 119 and 
(X 2 , Y 2 ) represent the coordinates of region 123. 

Thus, (X r Y x ) = (123°, - 0.062), 

(X 2 , Y 2 ) = (137°, - 0.188) 

and (X, Y) = (130°, Y), since AZI is the independent variable in this case. The 
deviation Y is computed as follows: 
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Table 12 


SECTORS 


Sector Range Window Geographic Azimuth 

(sec/Km) Window (degrees) 


A 

0.035 < 

U< 

0.052 

320 < 0 < 45 
(Through 360) 

B 

0.060 < 

U< 

0.070 

45 < 0< 75 

C 

0.067 < 

U < 

0.090 

110 < 0 < 170 

D 

0.036 < 

U C 

0.067 

140 < 0 < 170 

E 

0.036 < 

U < 

0.090 

170 < 0 < 190 

F 

0.036 < 

U < 

0.048 

225 < 0 < 250 

G1 

0.039 < 

U < 

0.090 

290 < 0 < 304 

G2 

0.039 < 

U < 

0.090 

304 < 0 < 312 

G3 

0.039 < 

U < 

0.090 

312 < 0 < 320 


Regions Contained 
in Sector 

112, 113, 114, 115 
116, 117 

118, 119, 120, 121 
122, 123, 124 

125, 126, 127, 128 
129, 130 

131, 132, 133 

134, 135 

103, 104, 110 

102, 105, 109 

101, 106, 107, 108 
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Figure 78, Relationship Between Regions and Sectors 
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Table 13 


LINEAR INTERPOLATION FORMULA 


Given: Points (X r and (X„, Y.,) 


The straight line Y = aX + b which passes through (X^, Y 
may be expressed as: 


Y = Y + 


X - x 1 

X 2~ X 1 


( Y 2 - Y l> 


j) and (X 2 , Y 2 ) 
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Table 14 


SUBARRAY FORMULAS FOR SECTORS 
SUBARRAY B1 


Sector 

Region 

U 

(Sec/Km) 

Azimuth 

(degrees) 

Deviation 

(seconds) 

A 

112 

0.0435 

348 

0.144 


113 

0.0435 

355 

0.147 


114 

0.0435 

22 

0.210 


115 

0.0435 

39 

0.147 

B 

116 

0.065 

52 

0.090 


117 

0.065 

73 

0.073 

C 

118 

0.078 

119 

-0.011 


119 

0.078 

123 

-0.062 


123 

0.078 

137 

-0.188 


122-124 

0.078 

147 

-0.209 


121 

0.078 

154 

-0.253 


120 

0.078 

161 

-0.273 

D 

125 

0.064 

145 

0.238 


126 

0.060 

145 

0.184 


127 

0.057 

145 

0.131 


128 

0.054 

145 

0.078 


129 

0.050 

145 

0.024 


130 

0.046 

145 

0.066 

E 

131 

0.085 

185 

-0.159 


132 

0.061 

185 

-0.075 


133 

0.054 

185 

-0.012 

F 

134 

0.042 

224 

0.151 


135 

0.042 

224 

0.151 

G1 

103 

0.077 

297 

0.054 


104 

0.070 

297 

0.098 


110 

0.043 

297 

0.044 

G2 

102 

0.082 

308 

0.041 


105 

0.067 

308 

0.031 


109 

0.048 

308 

0.016 

G3 

101 

0.081 

316 

0.100 


106 

0.064 

316 

0.001 


107 

0.058 

316 

-0.053 


108 

0.054 

316 

0.011 
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(X - x x ) 


Y Y 1 + (X 2 - X x ) ' (Y 2 " Y l* 


Y = -0.062 + (-0-188 - (-0.062)) = -0.125 seconds. 


III.2.6 Uncovered Areas 

It is also necessary to determine deviations for events which occur in areas 
not contained in any sector. Since such areas are not near areas in which numerous 
events are observed, little knowledge of the expected phase delays exists. One 
possible solution is to steer by means of plane waves with no subarray deviation 
correction. A second possible method is to interpolate or extrapolate using sets 
of regional averages from several regions. A third possible solution is to define 
more complex functions which will cover the entire teleseismic zone. This method 
appears most promising, but is contingent upon the results achieved by the pres¬ 
ently defined nine sectors. 


III.2.7 Testing The Regions And Sectors 

This section discusses a test which was performed to evaluate steering delay 
correction factors based on the construction of regions and sectors. The results 
of the test indicate that regional correction factors can be employed successfully 
to obtain improved phase delays. 

The test commences with a simulated steering of the array to an event located 
within a region, by means of plane wave delays modified by correction factors 
associated with the region. Next, a best fit plane wave is obtained for the event 
by using the time of arrival at the center seismometer of each subarray as input 
data for the Least Squares Wavefront Program. The difference between computed 
and actual arrival time at the center seismometer for each subarray is next cal¬ 
culated. Finally the Seismic Steering Delay Anomalies Program is used to evalu¬ 
ate the difference (if any) between the steering correction factor and the best fit 
plane wave deviation at each subarray. The Anomalies Program calculates from 
the 21 differences the loss (in dB) that can be attributed to using subarray steering 
correction factors. The magnitude of this loss measures the effectiveness of the 
set of steering corrections associated with the region containing the event. Since 
a region is characterized by its set of steering corrections, this loss also gauges 
the reliability of the region itself. If the loss is unacceptably large (where unac¬ 
ceptably large would have been previously specified), the conclusion is that one or 
more of the following factors is responsible. 

a. The region is not well-defined, i.e., the correction factors are 
ineffective for steering the array. 

b. The event was inaccurately timed at several subarrays. 


148 




c. A grossly inaccurate reading was made at one or more subarrays. 

d. Some other parameter, such as depth, must be taken into account. 

Factors b, c, and d can be eliminated by rechecking seismometer records and 
repeating the computations as required. If a recalculation continues to yield 
an unacceptably large loss, it must be concluded that the steering correction 
factors are of no benefit in reducing mis-steering loss. Such a test procedure 
will be more meaningful and reliable if there are several events per region which 
can be tested. 


III.2.8 Test Results 

The test conducted in this quarter was based on a set of 56 events, each with 
center seismometer arrival times, made available by Seismic Data Laboratory, 
Teledyne,Inc. 

The events were processed through the Least Squares Wave Front Program 
which calculated inverse phase speed and direction for each event. As shown in 
Figure 79, it was found that 33 events fell inside regions and 23 fell outside. 

The results of the test for the 33 events which fell inside regions are shown 
in Table 15. Columns 1 and 2 are self-explanatory. Column 3 indicates the num¬ 
ber of center seismometers that were used in the curve fitting process. The 
fewer the number used, che less reliable is the fit. Inspection of the loss figures 
given in Column 4 shows that 26 of the 33 events suffered a loss of 1 dB or less 
at 1 Hz. Thus it may be concluded that regional correction factors would have 
helped steer to these events successfully. As for the remaining 7, it is seen 
from Column 3 that for 4 of them, the number of sensors used was 10 or less out 
of 21. This small number of sensors implies that the wave fitting is not to be 
regarded as valid or reliable because the deviations from the best-fit wave front 
are so large that they cannot be compared with regional averages. Consequently 
these events have been disregarded. 

Interpolation procedures were not available for application to the 23 events 
falling outside the regions. However, some of these events lay rather close to 
regions and provided an opportunity to examine how far regional deviations could 
be extended. Accordingly, it was decided to simulate steering to eight of these 
events using for the simulation the correction factors for the regions closest 
to the events. These eight events are shown in Figure 80; the closest regions are 
shaded in the figure. The results of this attempt to stretch regional corrections 
are shown in Table 16. Columns 1 through 4 are similar to Columns 1 through 4 
of Table 16. Column 5 contains relevant comments. As can be seen from Table 16, 
five of the eight could have been steered to successfully using regional averages. 

The three events for which regional average steering would have been unsuccessful 
have been labeled by a check (\J) in Figure 80. 
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Figure 79. 23 Events Located Outside the 35 Regions 
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Table 15 


TESTING OF THE REGIONAL AVERAGES USING 
33 EVENTS LOCATED WITHIN REGIONS 


Region 

Number 

Event 

Identification 

Number 

Number 

of 

Sensors 

Loss 
(in dB) 

101 

4840 

6 

1.65 


4850 

13 

1.41 


4880 

20 

0.35 


5000 

20 

0.55 


5020 

18 

0.37 


5050 

18 

0.54 

103 

4830 

17 

0.33 


4980 

17 

0.60 


5100 

20 

0.57 

104 

4635 

20 

0.87 


4655 

20 

0.96 


4680 

20 

0.20 


5145 

19 

0.30 


5155 

18 

0.35 

106 

4625 

20 

0.51 

108 

4310 

19 

0.23 

110 

4639 

9 

1.0! 

114 

4930 

20 

0.70 

115 

4780 

16 

1.01 


4810 

14 

0.43 


4860 

17 

0.57 


4940 

21 

0.50 

116 

4627 

10 

1.41 

121 

4700 

11 

0.74 


4790 

19 

0.81 


4870 

10 

1.31 

124 

4920 

21 

0.26 


4950 

21 

0.32 
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Table 15 
(Concluded) 


Region 

Number 

Event 

Identification 

Number 

Number 

of 

Sensors 

Loss 
(in dB) 

125 

5060 

21 

0.37 

128 

4630 

21 

0.23 

131 

4800 

19 

0.34 


5010 

12 

1.45 

135 

4320 

20 

0.36 
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Figure 80„ Steering to 8 Events Near Regions 
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Table 1G 


TEST OF STEERING TO EVENTS LOCATED NEAR 



REGIONS USING REGIONAL AVERAGES 


Region 

Event 

Number 

Loss 

Comment 

Number 

Identification 

of 

(in dR) 



Number 

Sensors 



101 

4960 

21 

1.38 


102 

4623* 

11 

1.61 

This set 
unacceptable 

105 

4623* 

11 

0.71 

This steering 
acceptable 

109 

4629 

13 

0.23 



5115 

20 

0.23 


116 

4631 

17 

0.52 


118 

5080* 

18 

i. ig'I 

Neither set is 
acceptable 

119 

5080* 

18 

1.73 f 


129 

5070* 

18 

1.01 

Unacceptable 

130 

5070* 

18 

0.79 

Acceptable 

132 

4637 

15 

3.24 



*Event is steered to using two or more sets of regional averages. 
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III.2.9 Conclusions 


The test results verify that the present regions are well-defined. Correction 
factors associated with the regions have been used to steer successfully to nearlv 
all events located within the regions. The fact that five events located near 
regions were also steered successfully by using the regional corrections lends 
credence to the planned method of interpolating between regions. 
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III. 3 SUB ARRAY ASSIGNMENT 


III. 3.1 Introduction 

The task addressed here is to cover an arbitrary location, Uj , in inverse 
velocity space by a judicious choice of pre-formed, pre-steered subarray beams. 
Ordinarily the problem is resolved by choosing that set of beams, one from each 
subarray, which is closest to the point Uj. However, consider the events plotted 
in Figures 81, 82, and 83. The square ( □) in each figure denotes the location of 
the event in question. Observed seismometer arrival times were used to best-fit 
a plane wave to each of the 21 subarrays. ^Then for each subarray, the best-fitting 
plane wave was used to calculate a value U^j which may be interpreted as the 
event location as estimated by using only subarrayi. For each event, the 21 
values Ugj have also been plotted in Figures 81, 82, and 83_where they are shown 
as dots (.). It is seen from the figures that in no case is equal to U- for any 
pair (j*£). 

The situation depicted in these figures has several implications for beam¬ 
forming and subarray steering. Since the dots denote the positions to which the 
subarrays should be steered to form the best LASA beam covering the correspond¬ 
ing square (□ ), it is possible that for one or more subarrays the preformed beam 
which is closest to the nominal location Uj is not the same as the beam which is 
closest to U^j. If this is true, the assignment of subarray beams is no longer an 
elementary task. 

The fact that Uj and U^j may not be the same point indicates that in some 
cases phase delays will be obtained by a chaining process. The location being 
steered to will not supply the delays, but it will point to several other locations, 
each of which will furnish its associated delay. Organization of the phase delay 
library must take cognizance of this possibility. 


III.3.2 Notation and Definitions 

Notation for indexes,as used in this section, is defined as follows: 


NOTATION 

DEFINITION 

RANGE 

i 

instrument index 

1, • . • , I 

j 

LASA arr&y beam index 

1, • • • , J 

k 

subarray beam index 

1, . . . , K 

i 

subarray index 

1, .... E 
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Figure 81 „ Kamchatka Event 
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Figure 82„ Kazakh Event 
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Figure 83 o Novaya Event 
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The following bounds may be assumed for the above indexes: I £ 25, 400 £ J £ 1024, 
K £ 10, and L < 21. The inequalities I < 25 and L < 21 imply that selective assign¬ 
ment of instruments and subarrays may be utilized for beamforming processes. 


Predefined sets of points in inverse velocity space, assumed as inputs for 
the problem of subarray assignment, are as follows: 


a. The arbitrary set {Uj = (U x j, Uyj) } , assumed to be covered 

by a corresponding set of LASA beams. 


b. The set { U = (U Xj ^, Uy^j)} , where U^j is interpreted as an 
"aiming point" covered by that beam of subarray i which best 
contributes to the formation of LASA beam j, under optimal 
processing. This set of points is referenced to the set{ Uj} , 
each point being an estimated location of TTj, as viewed from 
subarray f. It is noted that in the absence of anomali 


subarray 
for each pair i, j. 


alies, U.: = TT.. 
J ij 


c. The set {Uf k = U i k )}, where Uf k is interpreted as an 

"aiming point" covered Dy the K^h beam of subarrayi. This set 
of points is assumed to have been selected on the basis of some 
previously defined optimization criteria and is subject to reassign¬ 
ment on the basis of these criteria. 


III.3.3 Evaluation of Loss 


Evaluating the maximum loss observed for any LASA beam is the first 
problem considered. This loss is associated with a specified steering pattern 
{Ui k } of preformed beams. A procedure for evaluating such loss is dis¬ 
cussed below and presented diagrammatic ally in the flow chart of Figure 84. 


First a specified LASA beam j is considered. For each subarray i the sub¬ 
array beam k is determined so that the vector magnitude At j k = - Uf k I i R 

minimal, and this minimal value is designated by the symbol 




The minimization is accomplished by calculating the set of numbers 
{^ijk (k = 1. 1 • • » K)} and selecting its smallest member A^*'. 


The program now calculates the subarray loss 
A k 0 , utilizing Formula 2.3 1 . 


i j 



= A^(f A^o ) 2 . 

ij 



corresponding to 


In the above formula, f is frequency and is a constant which depends on the 
subarray geometry. 
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Programming 
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ng order, where 

LL. = Max.{ LL.} 
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I 


Future 

Programming "1 


Calculate LAS A 
Lass 

LL. 


Add 1 ta i 


J 


Output 


A. List af LASA Beam Locations {IL } 

and Corresponding Subarray 
Aiming Paints{ (T . } 

AJ 

B. List af Subarray Beams {U^} Farmed 
far LASA Beam Assignment 

C. l. List af Assigned Subarray Losses 

and Corresponding Subarrays, Sub- 
arrays, Subarray Beams, and LASA 
Array Lasses 

2. List of Maximum Subarray Losses 
and Corresponding Subarrays, Sub- 
array Beams, and LASA Array Lasses 

3. Sorted List af LASA Losses {LLj } 

D. List af Subarray Beam Assignments 
and Assignment Frequencies 

Far Each U . Supply the Corresponding 
4ka * a Delay Generation 

Program 


. Yes 

Calculate 


ii 

£ 

5* 

7T 

7T 


Calculate Sub¬ 
array Lass 

\j * k v 


Add 1 ta i 


a y *> 

/»\ 

✓ LL. x 

: jo > 

Minimal / 

X s 


Na 


_1_ 

I I 

I f - , I 

j Reassign {U^} | 

I ' 


Fsgure 84. Subarray Assignment Program 
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There is a loss at the LASA level, represented by the symbol LLj, which 
occurs because of the use of pre-formed beams. This loss is a function of the 

L subarray losses: {X ... } (i = 1, . . . , L). Loss LL. is given by: 

J 


LL. = -20 log in [ T 



/ 20 > 


Subarrays will be identified by number so that if one or more are deleted, their 
numbers will be passed over in forming LLj. The purpose of such numbering 
is to permit identification of exactly which subarrays are being used in cases 
where the LASA beam is formed using fewer than the full number of subarrays. 

The preceding calculations are iterated to obtain the sets of subarray losses 
{ ^j£k 0 i } (j = 1, . . . , J; i = 1, • • . , L) and the set of LASA losses {LL.} 

(j = 1. ! . . , J) associated with all LASA beams previously defined. The •* 

LASA losses {LL- j will be output in decreasing order. In addition, the computed 
subarray losses ^ will be listed. An independent output of the above 

program will be a list specifying the subarray beams to be used in forming 

each LASA beam U».. This list will be supplied to the Time Delay Generation 
Program. ^ 


III.3.4 Minimization of Loss 

Remaining problems that must be formulated more precisely and for which 
solutions are needed are as follows; 

a. Formulation of criteria for acceptable minimization of the maximum 
loss LLj 0 of the set {LLj } (j = 1, . . . , J). 

b. Formulation of criteria for iterated reassignment of the set { Ujfo } of 
subarray assignments. Such criteria should provide a method of 
feedback control maximizing the probability that iteration of the 
outer loop of Figure 84 will produce successively smaller values 

of LL; . 

Jo 
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Appendix IV 


EVENT PROCESSING STUDIES 


Event processing, as outlined in Reference 4, uses seismometer data along 
with detections supplied by the Detection Processor to estimate the location, depth, 
magnitude, and origin time of an event. Inherent to the system will be the neces¬ 
sary filtering, recording, and display processes with optional operator intervention. 

This appendix contains two studies which support event processing. The first 
is concerned with an energy "density" method of calculating event magnitude. The 
method is related to the classical measure of magnitude for a sine wave ground 
displacement model. The second deals with the least squares orthogonal poly¬ 
nomial fitting method used to reduce voluminous discrete travel time data to a 
polynomial form. The polynomials can be used by the Event Processor to deter¬ 
mine post P arrival windows and to estimate geographic locations of an event. 


IV. 1 MAGNITUDE ESTIMATION 

An integral part of event processing will be the capability 7 to automatically 
estimate the parameter magnitude which is related to the energy released at the 
source of the event. The current classical technique uses manual measurement 
to obtain waveform amplitude and dominant period. Computer extraction of these 
parameters appears difficult especially when the signal-to-noise ratio is small. 

In this section an average energy estimation method is presented and is 
related to the calculation of event magnitude for a sinusoidal model of ground 
motion. Furthermore, the method is expected to produce reasonable estimates 
of magnitude for non-sinusoidal waveforms because, unlike the classical method, 
all spectral components contribute to the average energy estimate. 


IV. 1.2 Classical Determination of Magnitude 

18 

Gutenberg and Richter have related the unified "magnitude" of a seismic 
event to the compressional body wave amplitude received at a seismic station as: 

10^A 

m b = Q ( A > h ) + lo S 10 -Ip- + S, (89) 
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where Q(A,h) is the depth distance function for P wave using a vertical instru¬ 
ment, A is the maximum zero to peak amplitude of ground displacement in 
the first three cycles and measured in meters, T is the dominate period in 
seconds, S is the station correction, A is the epicentral distance in degrees 
between the event and the seismic station and h is the event depth in kilometers. 

The above formula gives consistent results for manual measurements of A 
and T from seismograms. It is, however, not easily implemented on the digital 
computer, nor does it account for the spectral characteristics of the signal and 
difficulty is expected when the signal to noise ratio is low. 


IV. 1.3 Sine Wave Ground Motion Model and Magnitude Estimate 
19 

Richter ' notes: "The quantity A/T has the advantage of being simply 
related, in theory at least, to the kinetic energy of its wave train." Using this, 
we assume that the ground displacement along the axis of a seismometer result¬ 
ing from a distant event is of the form: 

g(t) = A sin — t meters, (90) 


where A is the maximum ground displacement in meters, T is the period in 
seconds and t is the time in seconds. The ground velocity is: 


v(t) - 


2 7rA 

T 


cos 



meters/sec. 


(91) 


The kinetic energy of an incremental mass beneath the seismometer is: 

E’(t) = 1/2 6 m v^ (t) joules, (92) 

where 6 m is the incremental mass in kilograms. The average kinetic energy 
normalized with respect to incremental mass is: 

< E(t) > = TT 2 (y) 2 joules/kg. 

The average energy for sinusoidal ground displacement can now be related to 
event magnitude. 

To relate the average energy measurement of the sinusoidal model to 
magnitude, solve Equation (92) for | A j and substitute in Equation (89) to give: 

m b = 0.5 1og 1() | < E(t)> | + Q(A, h) + K + S, and (93) 

K = 5.502502 = -log Q 7r + 6.0 (94) 


164 




Equation 93 is valid for a sine wave input and should give a reasonable estimate 
of magnitude for other periodic waveforms because it includes the average power 
contribution of all spectral components. However, it does not include the effects 
of seismic noise. 

Because low level seismic noise is always present, the average energy 
measurement indicated in Equation (92) is a sum of signal and noise energy. 
Assuming that the noise is stationary and random: 

<E s+n (t)>= <E s (t)>+ <E n (t)> ’ ( 95 ) 


where E c ,(t) is the kinetic energy of the signal per unit mass and E n (t) is the 
kinetic energy of the noise per unit mass; so that 


<E (t)> = <E , (t)> 
s' ’ s+n' ' 


<E (t)> 
n' ' 


(96) 


Therefore, average signal energy can be found by subtracting the average noise 
energy from the measured average signal plus noise energy. 

The magnitude estimation process is outlined in Figure 85. The signal is 
not truly periodic, but rather it is pulse-like in nature and has a large amplitude 
for just a few seconds of time. It is during these few seconds that the signal plus 
noise average energy must be estimated. The magnitude estimation process is 
accomplished by sliding a window along the waveform and calculating the average 
energy in the window. When the window covers the signal, its output is maximum. 
The peak value is taken as the estimate of the signal plus noise average energy. 
Equation (93) is rewritten as 

m b = 0.5 1og 10 I <E g (t) > I p + Q( A, h) + K> S (97) 

where <E^(t)> is calculated from Equation (96) at the peak value of <E g+n (t)>. 

The signal plus noise energy estimation window should be the same length as the 
significant signal. This is about one or two cycles depending on the signature of 
the event. Similarly, Equation (92) cannot be evaluated exactly for seismic noise 
because it is stationary for only a few seconds. Therefore, the average signal 
energy and consequently the event magnitude, is based on the two estimates of 
the average noise energy and the average signal plus noise energy. 


IV. 1.4 Sensor and Transfer Function 

The sensor system transfer function is depicted in Figure 86. The velocity 
transfer function is much flatter than the displacement transfer function indicating 
that the seismometer is essentially a velocity measuring device. Using the instru¬ 
ment's velocity measuring capability, a direct approximation of the quantity A/T 
can be obtained and employed as indicated above. 
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Figure 85. Automatic Magnitude Estimation Process 
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Figure 86. Instrument Response 






IV.2 TRAVEL TIME AND HORIZONTAL INVERSE PHASE VELOCITY 
CHARACTERIZATION 

Seismic body wave travel time along with the related horizontal inverse 
phase velocity are featured in two areas of event processing.4 First, the 
estimation of post P arrivals in both time after the P arrival and also 
location in the horizontal inverse phase velocity space (hereafter designated 
as U space). The second area of interest is the conversion of U space loca¬ 
tions and arrival times for P and post P waves to geographic coordinates and 
an estimate of event depth. 

The purpose of this section is to demonstrate the use of least squares 
orthogonal curve fitting in one dimension to reduce the necessary travel time 
and U space magnitude data to a minimum. Curve fitting also appears to auto¬ 
matically assist in interpolation and in obtaining the derivative of the discrete 
data. 


IV.2.1 The Jeffreys-Bullen Travel Time Tables 

15 

The Seismological Tables by Jeffreys and Bullen , are used to determine 
event location and time of origin from event arrival times recorded at a set 
of seismic stations. They have been employed to aid in identifying a method 
for concise data representation and, quantitatively, how far such a process 
need be carried out. 

The total travel time, T, for a seismic body wave can be found directly in 
the J-B tables as a two dimensional function of event range and event depth. 

The magnitude of the horizontal inverse phase velocity vector in seconds per 
kilometer is by definition the partial derivative of the travel time with respect 
to distance. 

Finding each value of |u| from the J-B tables would require calculating 
the derivative of a low order smoothing function passed through several imme¬ 
diate data points or the physical storage of a corresponding table. Storage and 
manipulation of large tables (P alone has over 1400 entries) can be a burdensome 
process. Direct usage of the J-B tables to find total travel time, T, and the 
inverse horizontal phase velocity from an assumed range and depth is a case 
in point. While the tables are consistent in having 14 fixed depths, the epicentral 
range increments are not the same for all arrival waves. In particular, the 
change of range from one entry in the table to the next varies considerably in the 
depth allowance tables. 

As an aid in identifying later arrivals of an event, there is a need to compute 
theoretical travel time differences between an assumed P wave and the later 
arrivals, and the utilization of approximating analytical functions must be 
considered. 
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IV.2.2 Table Approximation by Classical Least Squares 

Since direct usage of the J-B tables would be computationally unwieldy in 
event processing, attention was given to the classical curve fitting of data points 
in the least squares sense. The fit to P wave data of surface depth deteriorated 
beyond fourth order despite the use of double precision (16 digit) arithmetic 
computer operations. This was not surprising, since the classical method 
requires the inversion of a matrix and involves the notorious principal minor of 
the well known Hilbert matrix^. 

Curve fitting, by the use of polynomials which are orthogonal under summa¬ 
tion, was deemed the next logical step since their use does not require the inver¬ 
sion of a matrix. While tabulations of precalculated orthogonal polynomials are 
available*^, their application requires equal spacing of the independent variable 
increments. 

Since the J-B tables have an unequally spaced variable of epicentral range, 
Forsythe's regression form of orthogonal polynomial curve fitting was adopted“^*>22. 
A brief description of this method follows. 

Given data for (m) points (x k , y k ), k = 1, 2, . . . , m; express (y) as a func¬ 
tion of (x) by the use of orthogonal polynomials of degree 0 through degree (n), 
where (n < m - 1). Thus, find bg, bj, b2» ..... b n so that: 

y(x) = t> 0 P 0 (x) + bjP^x) + . . . + b n P n (x), (98) 

where P^(x) is a polynomial of degree (i) and y(x) best fits the data in the sense 
of least squares: 



- y ( x k )l 


2 . 

is a minimum, 


where y k is the given tabular value and y(xk) is the corresponding calculated value. 


After y(x) has been expressed in the form 98, it can be converted to a poly¬ 
nomial in powers of x since each Pi(x) is a polynomial in x: 


2 n 

y(x) = a n + a x + a x + . . . + a x 


(99) 
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The coefficients b (( , b^, b.,, . . . , b n may be found from the formula: 


b. = 

1 



m 




k=i 


( 100 ) 


The orthogonal polynomials are generated by a recursive formula as 
follows: 


P^x) = 0 


P 0 (x) 


P. +1 (x) = X P.(x) - fc i+1 P.(x) - ^ Pj.^x). 


( 101 ) 


The coefficients and (3 . are chosen so that the polynomials P^, P^, P 
, P n are orthogonal under summation: 

m 

\ p i( x k ) p j( x k) = 0 for a11 ( i ^ J)- < 102 ) 

k-1 


The formulas for a. 1 and B. are as follows: 

l+l ' l 


m 


I x k (P i 


(Xu))‘ 


a. , = 


k=l 


1+1 m 


X*. 

k-1 


( x i,)) 


(103) 
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m 


( p i(x k >r 


p . = 


k=l 


i m 


( p i_l(x k ))' 


<=1 


(104) 


Forsythe suggests transforming the independent variable x to the range -2. 
to +2. to avoid possible overflow. This transformation also permits fewer 
digits to be used for each coefficient a^ for a specific accuracy of fit. 

IV.2.3 Results of Least Squares Orthogonal Polynomial Fit to the 
Jeffreys-Bullen Tables 

The Forsythe regression form of orthogonal polynominal curve fitting was 
used to reduce the tables to manageable polynomials. The errors of fitting the 
travel time, T, for the P wave, PcP wave, and pP-P wave difference arc evident 
on figures 87 to 90. 

Figure 87 is primarily a plot of the error of fit in seconds versus the order 
of the polynomial fit. Two indications of fitting error are given: the standard 
deviation (RMS error) and the maximum deviation. Both are plotted as a range 
covering the 14 depths. Using | U | - 0.04 sec/Km at 30° distance as a worst case, 
a timing error of one second is equivalent to a 25Km error in location and is 
called the ’location error bound” on the graph. Q is the quantization level of the 
original J-B data and a is the expected RMS error in the original data caused 
by quantization. 

The ninth order fit for T seems best because the maximum location error is 
only 3 Km and the RMS error value of the fit was bounded by a. 

Figure 88 is very similar to Figure 87 except that the wave in question is 
PcP and for an event at the surface only. Using |u| = 0.023 sec/Km at 30° dis¬ 
tance as a worst case, a timing error of one second is equivalent to a 43.5 Km 
error in location. An eighth order fit looks best for T. 

Figure 89 shows the error of fit for the 13 depth correction polynomials for 
PcP. The d/RQ is a standard form for showing depth in a compact form and can 
be viewed as ratio of the event depth below the crust (the crust is about 33 Km 
thick) to the radius of the earth minus the crustal layer. Thus, a depth of 
33 Km, as measured from the surface of the earth, becomes equal to zero in 
d/Rg format. The worst case for an error in the estimation of depth due to a 
timing error occurs at 80° range and d/Rp = 0.11. Here, only 5.1 seconds 
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Locotion Error Bound in 
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Location Error Bound in Km 
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Depth Error Bound in 






























separates the two depth columns (0.11 and 0.12) in the J-B tables. This results 
in an upper bound of 12.4 km error in depth for a one second timing error. 

Figure 90 is similar in design^to Figure 89 and shows the error of fit for 
the pP-P time difference. Since | U | is the derivative of travel time with respect 
to distance, the order of the |u | polynomial will be one order less than the proper 
fit for T, 


IV.2.4 Conclusions and Recommendations 

Forsythe's regression form of orthogonal polynomial curve fitting outper¬ 
forms the classical least squares method and appears to be quite useful in fitting 
one dimensional travel time data. While taking the derivative of the polynomial fit 
to obtain | 0] as a function of range seems better than the derivative of discrete 
data, it may offer only marginal accuracy in geographic conversion yet be per¬ 
fectly acceptable for post P detection windows in U space. Fitting T and | u| data 
obtained at LASA should be no problem since they are direct measurements. It is 
recommended that a two dimensional surface fitting method be developed so that 
both range and depth are independent parameters. 
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